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Abstract 

The  near  future  will  see  a  proliferation  of  small  low  cost  communication  and  science  satel¬ 
lites  with  modest  (~0. 1  -  0.5  deg)  pointing  requirements  which  will  use  attitude  determina¬ 
tion  systems  (ADS’s)  of  low  power,  weight,  size  and  cost.  This  thesis  is  a  study  into  the 
integration  issues  of  such  a  system,  that  uses  micro-mechanical  inertial  sensors  and  an 
interferometric  GPS  (IGPS)  attitude  determining  receiver  that  are  tightly  coupled  using 
“model  based”  estimation  techniques.  The  integration  issues  analyzed  include:  the  effect 
on  performance  of  limited  IGPS  use,  the  effect  on  performance  of  various  gyro  qualities 
and  model  accuracies,  the  inherent  robustness  to  failure  of  the  system  and  a  power  con¬ 
sumption  analysis  for  a  minimum  energy  IGPS  receiver. 

It  was  found  that  attitude  performance  varied  less  than  0.05  degrees  (RSS)  when  the  inter¬ 
val  between  GPS  phase  measurements  was  varied  from  once  every  second  to  once  every 
30  seconds.  This  justifies  the  development  of  a  minimum  energy  receiver  to  take  advan¬ 
tage  of  the  power  savings  of  less  frequent  GPS  measurements.  During  times  of  low 
dynamic  model  accuracy,  the  IGPS/MM  Gyro  ADS  was  bound  by  the  steady  state  attitude 
error  of  the  IGPS  which  averages  a  1  a  RSS  error  of  about  0.28  degrees  within  the  GPS 
update  range  mentioned  above.  Attitude  error  also  depends  on  the  gyro  quality,  and  vari¬ 
ous  projected  system  accuracies  are  given  using  several  gyro  error  models.  Real  micro¬ 
mechanical  gyro  data  from  a  prototype  gyro  was  successfully  integrated  with  real  GPS 
phase  measurements  using  an  Extended  Kalman  filter.  Failure  scenarios  and  power  expen¬ 
ditures  were  also  analyzed  using  steady  state  linear  covariance  analysis  as  well  as  a  space 
simulation  which  uses  orbital  parameters  of  a  telecommunication  satellite  constellation. 
Total  power  expenditure  of  200  mWatts  can  be  achieved  for  this  system  and  still  fulfill 
mission  requirements. 
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Chapter  1 


Introduction 

The  next  few  years  will  witness  an  explosion  in  the  demand  for  communication  infrastructure. 
With  the  rising  popularity  of  the  internet  and  an  ever  increasing  need  for  global  media  communi¬ 
cation  such  as  fax,  voice,  data  and  video  lines,  the  telecommunication  industry  will  face  the  for¬ 
midable  challenge  of  trying  to  meet  the  demand.  In  addition  to  fiber  optic  lines,  various  small 
satellite  constellations  are  being  proposed  and  designed  to  provide  large  numbers  of  small,  low 
cost  communication  platforms  in  low-Earth-orbit  (LEO)  and  medium-Earth-orbit  (MEO).  Some 
of  these  projects  include  Iridium,  Globalstar,  Teledesic  and  Ellipso. 

1.1  Objectives 

The  Charles  Stark  Draper  Laboratory  has  been  developing  techniques  for  a  new  attitude  determi¬ 
nation  system  (ADS)  for  the  small  satellite  market.  This  project  proposes  to  synthesize  two  in- 
house  navigation  systems  as  an  alternative  for  the  traditional  LEO/MEO  spacecraft  ADS  suite. 
This  system  consists  of  the  Draper  Micro-Mechanical  (MM)  gyroscopes  aided  by  an  interfero¬ 
metric  GPS  (IGPS)  receiver.  The  two  systems  are  integrated  in  a  tightly-coupled  design  through 
an  Extended  Kalman  filter  (EKF)  described  in  Chapter  3. 

Mechanical  wheel  gyros  have  been  extensively  used  on  spacecraft  ADS  in  the  past.  The  main  dis¬ 
advantages  of  these  gyros  are  their  relatively  large  weight,  power  and  volume  requirements;  their 
cost;  and  their  need  for  measurements  from  other  sensors  to  bound  error  growth  and  keep  their 
biases  calibrated.  The  size  and  cost  issue  can  now  be  overcome  through  advancements  in  micro¬ 
machining  and  Application  Specific  Integrated  Circuits  (ASICs),  which  have  led  to  the  develop¬ 
ment  of  the  MM  gyros.  The  error  growth  and  calibration  challenge  will  be  met  by  integrating  an 
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attitude  determining  IGPS  receiver  with  the  gyros.  The  IGPS/MM  Gyro  ADS  will  use  low  quality, 
high  output-rate  gyros  to  measure  dynamics  while  the  interferometer  will  serve  to  bound  errors 
through  measurements  at  lower  sampling  frequencies. 

As  with  any  design,  the  proposed  system  must  meet  cost,  performance  and  reliability  criteria.  This 
thesis  addresses  some  of  these  issues.  Increased  performance  usually  comes  at  a  price,  and  this 
study  will  address  such  trade-offs  associated  with  this  system’s  integration.  An  analysis  is  pre¬ 
sented  in  Chapter  4  which  summarizes  the  effects  of  using  varying  IGPS  measurement  updates, 
varying  gyro  qualities,  and  varying  model  confidences.  Chapter  5  presents  a  demonstration  of  the 
robustness  to  failure  inherent  with  this  ADS.  Chapter  7  addresses  the  power  issue  associated  with 
the  system. 

1.2  The  Sensors 

A  brief  description  of  the  theory  and  hardware  pertaining  to  the  two  sensors  used  in  this  applica¬ 
tion  is  presented  below. 

1.2.1  IGPS 

The  use  of  GPS  interferometry,  in  which  the  differential  phase  between  antennas  is  used  for  atti¬ 
tude  determination,  is  an  increasingly  significant  application  of  the  GPS  system.  Several  commer¬ 
cial  three-axis  attitude  determination  receivers  such  as  the  Trimble  Quadrex  and  TANS  Vector,  the 
Ashtech  3DF  and  the  Loral  Tensor,  have  been  developed  and  marketed  as  off-the-shelf  technol¬ 
ogy.  The  commercial  IGPS  receivers  use  a  least  squares  point  solution  when  computing  attitude. 
Recent  interferometry  research  has  focused  on  using  Kalman  filtering  of  these  attitude  solutions 
to  take  advantage  of  previous  attitude  knowledge  and  gain  improved  performance. 
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Theory 

In  IGPS  the  phase  of  the  LI  signal,  the  primary  wave  carrier  transmitted  by  the  GPS  system,  is 
used  to  estimate  the  orientation  of  the  body  frame  with  respect  to  a  known  reference  frame.  The 
carrier  phase  is  tracked  and  differenced  between  two  antennas  to  produce  a  measurement  of  rela¬ 
tive  range  to  the  GPS  satellite  vehicle  (SV).  The  following  figure  illustrates  this  principle. 


Arij,  the  relative  range,  is  the  projection  of  the  baseline  vector  bt  onto  the  line-of-sight  (LOS)  vec¬ 
tor  p j  for  the  GPS  SV  on  channel  j: 

Arij  =  bJ  ■  P j  =  |M  •  cos(0)  (1.1) 

0  is  the  aspect  angle  between  the  bt  and  py- .  The  primary  observable  for  GPS  interferometry  is 
differential  phase,  Acp(y,  which  is  the  difference  between  slave  antenna  carrier  phase,  cp ,y,  and 
master  antenna  carrier  phase,  <pM  j.  Acp,-;  differs  from  delta  range  by  the  integer  ambiguity  A kiJt 
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which  is  some  integer  number  of  wavelengths  of  the  GPS  signal.  The  relationship  between  mea¬ 
sured  phase  and  differential  range  is: 


A  ru  =  A<Pu+Mi,rhi  d-2) 

where  etJ-  are  the  errors  in  the  received  signal.  These  errors  will  be  further  discussed  in  Chapter  2. 
From  equation  1.1  and  1.2  we  can  solve  for  the  measurement  equation: 


A%j  =  b,  Pj-A klJ  +  eu 


(1.3) 


The  method  for  determining  Aktj  used  in  this  thesis  is  the  decoupled  and  constrained  ambiguity 
resolution  method  (DCAR)  developed  by  Puri  [20],  Additional  information  about  IGPS  theory  is 
presented  in  Appendix  C. 


Hardware 

The  receiver,  which  collected  all  differential  phase  data  used  in  the  roof-top  filter  to  be  described 
in  Chapter  4,  is  the  Trimble  TANS  Vector  GPS  receiver.  The  Vector  is  a  six-channel,  four-antenna, 
C/A  code  receiver  which  outputs  differential  phase  measurements  at  the  rate  of  1  Hz.  The  six 
channel  tracking  loops  are  time  multiplexed  across  the  four  antennas.  The  receiver  is  connected  to 
the  four  micro-strip  patch  antennas  by  coaxial  cables.  The  patch  antennas  are  mounted  on  Trim¬ 
ble-supplied  8-inch  radius  metallic  ground  planes  to  reduce  susceptibility  to  multipath.  Each 
ground  plane  is  mounted  on  a  2  meter  maximum  length  adjustable  antenna  array.  This  array  was 
set  at  a  0.6  m  per  side  configuration,  yielding  a  total  of  three  baselines  from  the  master  antenna: 
two  60  cm  baselines  and  one  84  cm  baseline  across  the  diagonal.  The  array  was  mounted  on  a  2- 
axis  turntable  on  the  roof  of  the  Draper  Lab  building.  The  phase  data  was  collected  and  stored 
using  a  laptop  computer  which  was  connected  to  the  receiver.  The  following  picture  illustrates  the 
IGPS  hardware. 
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Figure  1.2:  IGPS  Data  Collection  Hardware 
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1.2.2  MM  Gyro 

The  development  of  micro-mechanical  tuning  fork  gyros  is  the  subject  of  current  research  and 
development  at  the  C.S.  Draper  Laboratory.  The  low  power  consumption  and  compact  size  of 
these  gyros  make  them  ideal  for  satellite  application. 

Theory  and  Hardware 

The  basic  operational  concept  of  the  MM  tuning  fork  gyro  is  the  sensing  of  a  coriolis  force.  Two 
proof  masses  vibrate  in  opposing  motion  perpendicular  to  the  axis  the  gyro  is  trying  to  measure. 
The  Coriolis  effect  is  described  by 

a  =  2m  xV  (1.4) 

where  a  is  the  acceleration  of  mass  with  respect  to  the  inertial  frame,  m  is  the  angular  rate  of  the 
gyro  frame  and  V  is  the  linear  velocity  of  mass  with  respect  to  the  gyro  frame.  The  following  fig¬ 
ure  shows  the  general  description  of  the  sensor: 
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Tuning  Fork  Gyro 


SIDE  VIEW 


w/  Motor  Motion 


Figure  1.3:  Illustration  of  MM  Gyro  Dynamics 

The  tuning  fork  gyro  consists  of  a  silicon  structure  suspended  above  a  glass  substrate  containing 
metallization  deposited  for  sensor  interfacing.  It  is  this  silicon  structure  that  suspends  the  two 
masses  by  a  sequence  of  beams  anchored  at  specific  points  in  the  substrate.  The  motor-drives 
force  the  two  masses  in  lateral,  in-plane  oscillatory  motion.  An  angular  rate,  go,  applied  about  the 
input  axis,  perpendicular  to  the  velocity  vector  of  the  masses,  generates  a  Coriolis  force  that  acts 
to  push  the  masses  in  and  out  of  the  plane  of  oscillation.  Since  the  velocity  vectors  of  the  masses 
are  in  the  plane  of,  and  orthogonal  to,  the  input  axis;  perpendicular  motion  orthogonal  to  the  plane 
is  induced  in  response  to  the  Coriolis  force.  The  resultant  motion  is  measured  by  capacitor  plates 
under  each  of  the  two  masses,  providing  a  signal  proportional  to  the  rate  input.  The  gyro  element 


25 


itself  is  very  small,  measuring  about  2  mm  on  a  side.  A  magnified  picture  of  the  MM  Gyro  is 
shown  in  Figure  1.4. 


Figure  1.4:  Micro-Mechanical  Tuning  Fork  Gyro 

The  bulk  of  the  size  and  complexity  of  the  MM  gyros  actually  comes  from  the  electronics,  which 
are  an  integral  part  of  the  gyro  system.  The  electronics  implement  the  control  loops  that  drive  the 
gyro  and  measure  its  out  of  plane  motion.  All  the  electronics  are  done  in  Application  Specific 
Integrated  Circuits  (ASICs)  using  Complementary  Metal-Oxide  Semiconductor  (CMOS)  pro¬ 
cesses  which  allow  for  the  entire  gyro  package  to  be  approximately  1  cm  X  1  cm  in  area  and  a  few 
millimeters  thick  [12] 


Figure  1.5  is  a  magnified  picture  of  the  MM  Gyro  and  ASICs  assembly.  It  measures  approxi¬ 
mately  one  centimeter  on  each  side. 
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Figure  1.5:  MM  Gyro  and  ASICs 


Chapter  2 


Error  Modeling 

2.1  Introduction 

The  performance  of  the  Extended  Kalman  Filter  that  combines  MM  gyro  measurements 
with  the  IGPS  measurements  inherently  depends  on  the  correct  modeling  of  the  system 
errors.  This  chapter  presents  models  of  the  IGPS  and  MM  gyro  measurement  errors  that 
will  be  used  in  the  Kalman  filter  design.  These  models  are  a  continuation  of  the  work 
already  done  in  IGPS/MM  Gyro  attitude  determination  at  the  Charles  Stark  Draper  Labo¬ 
ratory.  The  first  section  of  this  chapter  provides  models  for  the  error  sources  in  the  GPS 
differential  phase  measurements.  The  second  section  gives  a  standard  model  for  the  error 
sources  inherent  in  a  strapdown  gyro. 

2.2  GPS  Differential  Phase  Error 

This  section  summarizes  the  models  used  to  characterize  the  errors  present  in  GPS  differ¬ 
ential  phase  measurements.  Although  many  models  have  been  presented  in  the  literature, 
the  models  used  here  are  consistent  with  those  developed  by  Puri  [20]  at  Draper  Labora¬ 
tory.  The  emphasis  of  this  study  is  integration  issues;  hence,  the  models  already  developed 
were  deemed  sufficient  for  the  analysis. 

2.2.1  Overview 

Errors  in  the  GPS  differential  phase  measurements  can  directly  affect  attitude  accuracy.  A 
brief  description  of  the  sources  of  these  errors  follows.  The  mathematical  models  will  be 
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presented  in  the  next  section. 


Line  Bias 

The  LI  carrier  signal  arrives  at  each  patch  antenna  and  travels  to  the  Receiver  Processing 
Unit  (RPU)  along  four  separate  coaxial  cables.  Line  bias  error,  denoted  by  A(3,  is  any 
phase  delay  caused  by  cable  path  length  variations  from  the  antenna  to  the  receiver.  The 
bias  caused  by  this  delay  can  vary  since  the  path  lengths  may  change  with  temperature 
variations  of  the  system.  The  coaxial  cables  are  designed  to  be  as  close  to  the  same  length 
as  possible.  Assuming  the  cables  are  positioned  close  to  each  other,  any  thermal  variations 
which  affect  one  cable  should  affect  the  others  in  a  similar  fashion  since  attitude  is  deter¬ 
mined  by  phase  differences.  Thus  close  proximity  of  cables,  one  to  another,  should  help  to 
minimize  line  bias  error. 

Baseline  Length  Error 

The  attitude  solution  accuracy  is  limited  by  the  knowledge  of  the  baseline  vectors  (from 
master  to  slave  antennas)  that  define  the  vehicle  orientation.  The  baseline  vector  is  subject 
to  vehicle  flexure,  expansion  and  contraction.  On  a  spacecraft,  for  example,  one  source  of 
flexure  would  be  thermal  stresses  experienced  when  passing  in  and  out  of  eclipse.  Base¬ 
line  expansions  are  kept  small  by  placing  the  IGPS  array  on  a  rigid  platform  no  more  than 
one  meter  apart.  This  error  is  modeled  in  this  study  only  for  the  orbital  simulations  and  not 
for  analysis  performed  on  the  roof-top.  Roof-top  temperature  variations  were  not  great 
enough  justify  such  modeling.  Baseline  length  error  will  be  denoted  by  8b. 

Multipath 

One  of  the  largest  and  most  troublesome  causes  of  error  in  GPS  interferometry  is  multi- 
path.  Multipath  is  any  undesired  signal  reflection  that  interferes  with  the  direct  signal. 
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Since  the  reflected  signal  travels  a  longer  path,  its  phase  is  shifted  from  that  of  the  direct 
signal.  Ways  of  mitigating  multipath  include  limiting  any  reflectors  in  the  vicinity  of  the 
antennas  and  using  ground  planes  under  each  patch  antenna.  Once  multipath  has  already 
corrupted  the  signal,  there  are  ways  to  compensate  for  the  repeatable  part  of  the  error 
caused  by  static  sources  such  as  body-fixed  reflectors.  Cohen  and  Parkinson  [6]  have  dem¬ 
onstrated  cancellation  of  low  frequency  multipath  by  developing  a  spherical  harmonic 
model  from  experimental  data.  A  similar  fifth  order  model,  developed  by  Puri,  to  match 
the  Draper  roof-top  environment  will  be  used  in  this  study.  The  multipath  dynamics  not 
captured  in  this  bias  map  will  be  modeled  by  augmenting  the  receiver  correlated  noise 
model.  Multipath  is  not  expected  to  be  as  great  a  problem  on  the  spacecraft  as  on  the 
Draper  roof.  Total  multipath  error  will  be  denoted  by  Am. 

Antenna  Phase  Center  Variations 

Asymmetry  of  the  dielectric  substrate  of  the  antenna  which  captures  the  LI  signal  can  lead 
to  apparent  movement  of  the  phase  center  [24].  This  variation  is  a  function  of  the  inci¬ 
dence  angle  of  the  incoming  signal.  Since  typical  RMS  variation  of  2-4  mm  can  be  caused 
by  phase  center  variation  on  microstrip  GPS  antennas  [24],  it  is  important  to  calibrate  and J 
or  model  this  error.  Phase  center  variation,  A p,  is  body  fixed  and  time  invariant,  however,  it 
is  antenna  specific.  These  characteristics  allow  it  to  be  modeled,  along  with  low  frequency 
multipath,  in  the  calibration  model  mentioned  above.  Any  mismodeling  of  the  phase  cen¬ 
ter  by  the  bias  map  will  be  accounted  for  by  augmenting  the  receiver  correlated  noise 
model,  in  accordance  with  analysis  done  previously  [20], 

Receiver  Noise 

The  incoming  carrier  phase  is  corrupted  by  errors  induced  by  the  receiver  itself.  In  this 
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analysis  the  receiver  noise  is  broken  down  into  two  constituent  parts:  receiver  white  noise, 
Av;  and  receiver  correlated  noise,  A^.  The  white  noise  model  captures  all  of  the  uncorre¬ 
lated  errors  such  as  thermal  noise  and  oscillator  noise  [17].  Receiver  correlated  noise  mod¬ 
els  any  temporally  correlated  errors  in  the  receiver.  It  is  also  used  to  absorb  mismodeling 
in  the  multipath  and  phase  center  error  calibration  [20]. 


The  following  sections  present  models  for  line  bias,  baseline  length  error,  multipath,  phase 
center  variation,  and  receiver  noise.  They  are  all  in  accordance  with  [20]. 


2.2.2  Line  Bias  Error  Model 


The  single  difference  operator  subtracts  the  line  bias  of  each  of  the  slave  antennas  from  the 
master  antenna  bias  to  yield  the  differential  line  bias,  A(3: 


Ap, 

1-10  0 

AP  = 

a|J2 

ii 

10-10 

Ap3 

1  0  0  -1_ 

Pm 

P, 

P2 

P3_ 


(2.1) 


Since  the  line  bias  is  dependent  on  the  coaxial  cables,  the  differential  line  bias  is  a  function 
of  the  baseline  and  not  the  channel  from  which  the  measurement  comes.  Hence,  only  three 
line  bias  error  states  will  be  needed.  A  first  order  exponentially  correlated  error  model  is 
used  to  model  line  bias: 

-  Ap_0)  +  "-^(0  (2.2) 

where  T^p  is  the  time  constant  and  w^p  is  zero  mean  white  Gaussian  process  noise.  The 
nominal  line  bias,  Af30 ,  is  calibrated  before  estimation  [20].  The  time  constant  depends  on 
the  receiver  environment  and  is  modeled  as  half  the  orbital  period  for  the  orbital  simula- 
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tion  that  will  be  discussed  later  and  12  hours  for  the  roof-top  configuration.  The  driving 
noise  strength  is  a  function  of  the  steady  state  line  bias  covariance,  aAp, 


2  _  2'gAP 


lAP 


(2.3) 


and  since  it  is  a  product  of  the  single  difference  operator,  it  is  correlated  between  antennas. 
The  correlations  are  added  as  shown  below  for  a  single  channel  and  three  baselines: 


1  0.5  0.5 
0.5  1  0.5 
0.5  0.5  1 


(2.4) 


where  QAp  is  the  process  noise  covariance.  The  line  bias  steady  state  covariance,  aAp,  is 
the  value  that  needs  to  be  determined  and  the  method  for  so  doing  is  described  in  detail  in 
reference  [20]. 


2.2.3  Baseline  Length  Error  Modeling 

It  is  assumed  that  the  flexure  of  the  IGPS  array  will  be  limited  to  expansions  and  contrac¬ 
tions  of  the  baselines  due  to  thermal  variations  of  the  environment.  Another  first  order 
Markov  process  will  be  used  to  model  the  baseline  error  dynamics: 

=  -±-bb  +  wb(t)  (2.5) 


with  driving  noise  strength: 


2  2  ab 

=  T~ 


(2.6) 


The  baseline  length  error,  bb,  enters  the  differential  phase  as  a  function  of  geometry  and 
the  resulting  differential  phase  error  due  to  bb  is  presented  here  without  derivation  [20]: 
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(2.7) 


\j  =  <“£ '  fijWl 

where  /  is  the  baseline  and  7  is  the  channel  as  explained  in  chapter  1,  and  ub  is  the  unit  vec¬ 
tor  for  a  nominal  baseline. 

2.2.4  Multipath  and  Phase  Center  Error  Modeling 

As  mentioned  earlier,  the  deterministic  errors  produced  by  body  fixed  sources  and  phase 
center  variation  are  basically  time  invariant  and  only  a  function  of  the  GPS  line  of  sight 
vector  in  the  body  frame.  These  errors  can  be  compensated  for  by  comparing  the  measured 
phase  to  the  predicted  phase  for  a  static  receiver.  Puri  developed  a  bias  map  similar  to  the 
spherical  harmonic  model  developed  by  Cohen  in  which  the  sky  was  divided  into  parti¬ 
tions  that  individually  contained  corrections  for  that  line  of  sight.  However,  only  body- 
fixed  multipath  due  to  diffraction  off  the  ground  planes  could  be  calibrated  since  the  data 
collections  required  movement  of  the  IGPS  antenna  array.  Multipath  sources  on  the  roof 
of  Draper  Lab  include  cooling  stacks,  an  aluminum  shed,  and  a  skylight  from  another 
building.  The  multipath  might  be  a  result  of  planar  reflectors  or  multiple  point  reflectors. 
Any  uncorrected  multipath  or  mismodeled  phase  center  variation  is  modeled  by  increasing 
the  steady  state  covariance  of  the  process  noise  in  the  receiver  correlated  error  state  model. 

2.2.5  Receiver  Noise 

As  stated  earlier,  the  receiver  noise  has  been  separated  into  receiver  white  (wide  band) 

noise  and  receiver  correlated  noise.  For  this  study  the  differential  phase  wide  band  noise 
2 

strength,  cAv ,  is  experimentally  characterized  as  a  function  of  SV  line  of  sight  elevation  in 
the  body  frame.  The  receiver  measurement  noise  is  also  cross  correlated  via  the  single  dif¬ 
ference  operator: 
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where  0p  is  the  elevation  of  the  SV  on  channel  1  in  the  body  frame,  and  the  example  is 
for  two  baselines  and  two  satellites. 


The  time  correlated  errors  in  the  receiver,  as  well  as  the  multipath  and  phase  center  varia¬ 
tion  mismodeling  mentioned  above,  are  modeled  with  the  receiver  correlated  noise  param¬ 
eter.  A  first  order  Markov  process  is  used  to  model  this  error, 

=  -  ^"A§(0  +  h>a  §(0  (2.9) 


where  is  the  time  constant  and  is  white  Gaussian  process  noise.  The  driving  noise 
strength  is: 


2  _  2-qA,(ep/)2 


(2.10) 


where  the  steady  state  covariance  of  A^,  G^ ,  is  a  function  of  SV  elevation  in  the  body 
frame,  0p^ .  The  process  noise  strength  is  again  correlated  across  the  baselines,  as  shown  in 
this  example  for  a  single  channel  and  three  baselines: 


Qa^  = 


1  0.5  0.5 
0.5  1  0.5 
0.5  0.5  1 


2.2.6  Differential  Phase  Sensitivities 


(2.11) 


The  following  equation  expands  the  error  term  in  equation  1.3  to  include  the  above  mod- 
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eled  errors: 


A%  j  =  bi0 '  P j  ~  Aki,  j  +  AP«  +  Ebu  +  Ami,  j  +  APi,  j  +  Av«,  j  +  A^,  j 


(2.12) 


With  this  equation  we  can  determine  the  sensitivity  of  the  differential  phase  measurement 
to  each  modeled  error  source.  The  sensitivities  will  be  important  when  deriving  the  Kal¬ 
man  filter  measurement  sensitivity  matrix  in  the  next  chapter.  All  the  sensitivities  are 
straight  forward  except  for  the  baseline  length  error  which  is  obtained  by  taking  the  deriv¬ 
ative  of  equation  2.7.  In  the  order  discussed  they  are: 


dAcp,.  . 

3aP;  =  8(a) 

(2.13) 

^A(Pi,  j  -  T 

9(5 bk)  8( l,k)  p'-) 

(2.14) 
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3a  :i  - 

(2.15) 

9A  cp. 
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(2.16) 

9Acp,-  . 

3Av;;;  =  Kuwti) 

(2.17) 

s^(')  ■ 

(2.18) 

(2.19) 

2.3  MM  Gyro  Errors 

The  IGPS/MM  Gyro  ADS  treats  gyro  system  output  as  measurement  inputs  to  the  EKF  to 
determine  attitude  and  angular  rate  rather  than  as  an  instrument  that  determines  body  rate, 
as  is  normally  done  with  ADS’s  with  high  quality  gyros.  These  inputs  are  essential  for  the 
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model  based  ‘tightly-coupled’  filter  design,  alluded  to  earlier,  and  which  will  be  formally 
explained  in  the  next  chapter.  This  type  of  integration  is  possible  due  to  modeling  of  the 
dynamic  forces  such  as  the  gravity  gradient  torque. 

2.3.1  Overview 

There  are  several  sources  of  error  that  are  inherent  in  all  gyros  and  some  that  are  particular 
to  the  MM  gyros  currently  being  developed  at  the  Draper  Lab.  The  typical  sources  of  error 
are  captured  in  the  classical  gyro  model  presented  below.  Other  sources  such  as  scale  fac¬ 
tor  bias  and  misalignment  error  will  be  mentioned  as  will  errors  specific  to  the  MM  gyro. 

2.3.2  Classic  Gyro  Error  Model 

The  classic  gyro  error  model  is: 

SM(0  =  fflr(0  +  bg  (/)  +  vg(t)  (2.20) 

where  coM  is  the  measured  angular  rate,  CdT  is  the  actual  rate,  bg  is  a  time  varying  bias, 
and  vg  is  the  white  measurement  noise. 

Angle  Random  Walk 

The  white  noise,  v g,  is  also  known  as  drift  rate  noise  or  random  drift  rate  noise.  It  is  mod¬ 
eled  with  strength  Rg.  It  produces  “angular  random  walk”  (ARW),  created  by  integrating 
the  angular  velocty,  ,  over  a  short  interval  At.  The  variance  of  the  error  due  to  ARW  is 
given  by: 

°ARW  =  Rg'At  (2.21) 

where  At  is  the  integration  step  size  or  the  inverse  of  the  output  rate.  ARW  get  its  name 
because  the  noise  term  integrates  into  “random  walk”  in  the  angle. 
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Gyro  Bias 


The  gyro  bias  is  a  varying  drift  rate  in  the  gyro  that  can  cause  considerable  error  when 
using  low  accuracy  gyros  for  extended  periods  of  time.  One  benefit  of  the  IGPS/MM  Gyro 
ADS  is  that  the  differential  GPS  phase  measurements  will  constitute  a  way  of  periodically 
recalibrating  this  drift  rate.  The  gyro  bias,  with  units  of  angular  rate,  can  be  modeled  either 
as  a  Markov  process: 

|<v  =  -f(>»)+v*  <2-22> 

& 

or  by  a  simpler  random  walk  model: 

§fbs)  =  v„  (2.23) 


This  study  will  use  the  model  presented  in  equation  2.22.  The  strength  of  the  bias  noise  is 
modeled  by: 


(2.24) 


r\ 

where  a  b  is  the  variance  of  the  gyro  bias  steady-state  value  and  tg  is  the  time  constant 
associated  with  the  bias  dynamics. 


2.3.3  Other  Gyro  Errors 

In  addition  to  the  standard  angle  random  walk  and  gyro  bias  errors,  there  are  some  other 
errors  that  can  affect  the  quality  of  gyro  output. 

Scale  Factor  Error 

Scale  factor  error  arises  when  the  calibration  of  the  relationship  between  actual  angular 
rate  and  gyro  output  is  incorrect.  The  attitude  error  this  causes  is  proportional  to  the  prod- 
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uct  of  the  net  rotation  angle  with  the  scale  factor  error  [21].  When  real  gyro  data  is  used  in 
this  study,  it  will  be  assumed  that  scale  factor  error  has  either  been  calibrated  out  or  is  neg¬ 
ligible. 

Misalignment  Error 

Since  strapdown  gyro  packages  are  literally  strapped  onto  the  body  of  the  vehicle,  a  mis¬ 
alignment  error  can  result  if  the  transformation  relating  the  gyro  axes  to  body  axes  is  not 
known  perfectly.  For  this  application,  it  will  be  assumed  that  the  misalignment  errors  are 
either  calibrated  a  priori  or  are  small  enough  to  be  negligible.  Body  flexure  is  also 
assumed  negligible. 

Temperature  Bias 

The  MM  Gyros  currently  under  development  at  the  Draper  Laboratory  are  very  sensitive 
to  temperature  changes.  In  a  space  application  vehicle  electronics  could  easily  see  temper¬ 
ature  fluctuations  from  0°C  to  40°C  [14].  For  this  reason,  the  MM  Gyro  output  must  be 
temperature  compensated.  The  temperature  dependence  of  the  MM  Gyros  is  shown  in  Fig¬ 
ure  2.1  for  a  15  hour  drift  test: 
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Gyro  Rate  Measurement  Data 


Temperature 


Figure  2.1:  Raw  MM  Gyro  Data 

After  applying  temperature  compensation,  the  resultant  data  is  well  characterized  by  the 
noise-induced  bias  model  of  Equation  2.2: 


Residual  Gyro  Output  after  2nd  Order  Temperature  Fit 
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Figure  2.2:  Temperature  Compensated  Gyro  Data 


When  not  using  real  data  in  this  study,  gyro  measurements  were  created  by  simulation 
using  the  model  of  Equation  2.2. 
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Chapter  3 


Filter  Design 

3.1  Introduction 

The  key  to  the  tightly  coupled  integration  of  the  IGPS/MM  Gyro  ADS  is  the  Kalman  filter. 
It  combines  gyro  and  GPS  measurements  to  form  an  accurate  attitude  estimate.  As  men¬ 
tioned  in  Chapter  1,  this  application  requires  an  Extended  Kalman  filter  (EKF)  to  deal 
with  nonlinear  dynamics  and  nonlinear  measurement  relationships.  The  nonlinear  mea¬ 
surement  relationships  are  a  result  of  the  motion  of  the  GPS  satellites  that  provide  the 
phase  measurements. 

This  chapter  reviews  the  standard  filter  equations.  The  dynamics  used  in  state  propagation 
will  be  presented  along  with  a  formulation  of  the  measurement  relationships.  The  last  part 
of  the  chapter  will  give  a  brief  description  of  some  of  the  implementation  issues  which 
must  be  addressed  in  the  filter  algorithms. 

3.2  Filter  Equations 

A  review  of  the  Extended  Kalman  filter  equations  will  be  presented  in  the  following  sec¬ 
tions.  The  quaternion  formulation  of  the  state,  presented  in  Section  3.2.2,  was  first  demon¬ 
strated  by  Lefferts,  Markley,  and  Shuster  [15]  and  has  been  successfully  implemented  by 
many  researchers  including  some  in  the  IGPS  field  [3,  20,  25].  Additional  information  on 
Kalman  filtering  theory  can  be  obtained  in  Gelb  [9]  and  [1]. 
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3.2.1  Continuous-Discrete  EKF 


This  application  lends  itself  to  a  particular  class  of  estimation  problems  in  which  the  non¬ 
linear  system  has  continuous  dynamics  and  discrete-time  measurements.  The  system 
dynamics  are  characterized  by  the  continuous  nonlinear  state  differential  equation: 

x(t)  =  f(x(t),t)  +  g(x(t),t)w(t)  (3.1) 

and  measurements  are  characterized  by  the  discrete-time  nonlinear  measurement  equa¬ 
tion: 


h  =  h(xk’h)  +  vk  (3-2) 

where  the  subscript  k  refers  to  the  measurement  at  time  t*.  w(f)  is  a  zero  mean  white  Gaus¬ 
sian  process  noise  with  strength  Q(t): 

E[w(t)w(t-x)r]  =  Q(r)6(r-x)  (3.3) 

where  6  is  the  Dirac  delta  function  with  units  of  1/time.  vk  is  discrete  white  Gaussian  mea¬ 
surement  noise  with  covariance  R^: 


E[vkv]} 


Rjt  i=  k 
K  0  i*k 


(3.4) 


and  it  is  assumed  that  there  is  no  cross  correlation  between  the  measurement  noise  and  the 
process  noise: 


E[w{t)vTk]  =  0  for  all  t  and  k. 

We  define  the  error  between  the  true  state,  x,  and  the  state  estimate,  x  : 

X  =  x-x 

The  error  covariance  matrix  is  defined  as: 


(3.5) 


(3.6) 


P  =  £[(i)(i)7’]  (3.7) 

If  the  functions  of/  and  g  are  sufficiently  well  behaved,  we  can  linearize  the  system  about 
the  state  estimate  at  time  t : 
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(3.8) 


G(')-^sWO,  (3.9) 

yielding  the  linearized  dynamics, 

x(t)  =  ¥(t)x(t)  +  G(t)w(t)  (3.10) 

The  minus  and  plus  subscripts  denote  the  state  estimates  before  and  after  measurement 
updates,  respectively.  A  predicted  measurement  at  time  tk  can  be  formulated  by  using  the 
state  estimate  and  the  nonlinear  measurement  equation: 

zk  =  h(xk_,tk )  (3.11) 

Now  the  measurement  residual  is: 


zk  =  Zk-Zk  (3.12) 

which  can  be  related  to  the  state  error  vector,  x ,  by  the  following  linearized  sensitivity 
matrix,  Hfc: 

h  =  hA  +  v*  (3.13) 

where  H*  is  defined  by, 


H  =  —k 
k  dx 


(3.14) 


We  can  incorporate  the  measurements  to  form  the  optimal  state  correction  as  a  function  of 
the  residual  and  the  Kalman  gain,  Kk: 


xk+  =  (xk_  +  KkZk)  (3.15) 

and  the  updated  covariance  matrix  shown  in  the  Joseph  form: 

Pt+  =  (l-KtHt)P^  (l-K4H/  +  KtRtK[  (3.16) 

This  above  form  is  used  to  ensure  symmetry,  where 
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(3.17) 


Pt+  =  £[(*.„  )(it+  )r] 

and 

K^.=  Pk_  HJt(HA.Pjt_  H^+R^,)  (3.18) 

After  a  measurement  update,  the  state  is  propagated  to  the  next  time  increment  tk+] 
through  the  nonlinear  dynamics: 

x(t)  =  t )  (3.19) 

The  error  covariance  matrix  is  propagated  via  the  discrete  solution  of  the  continuous  time 
matrix  Ricatti  equation: 

P  =  FP  +  PFr  +  GQGr  (3.20) 

which  has  the  solution, 

P(*  +  D-  =  +  Pk+  <u  +  Q(^+u)  (3.21) 

where  O  is  the  state  transition  matrix,  and  Q  is  the  process  noise  dynamics  both  defined 
below: 


h  + I. 


+  =  j  F(x)0(t,  tk)dx 

lk 

(3.22) 

Qih’h)  =  1 

(3.23) 

h  + 1 

Q(tk  +  Vtk)=  j  <t>xG(%)Q(x)G(T)T<t>xdx 

1 1 

(3.24) 

where  1  is  an  identity  matrix  of  appropriate  size.  For  short  time  steps  these  equations  can 

be  approximated  to  first  order  by: 

®(tk+l,tk)  =  (1  +  FA  0 

(3.25) 

Q(**+i>f*)  =  G[QAt]GT 

(3.26) 
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where 


At  =  h  +  i-h  (3.27) 

3.2.2  Choosing  the  State 

A  critical  step  in  designing  a  Kalman  filter  is  deciding  what  states  to  incorporate  in  the 
system  dynamics.  The  states  of  interest  need  to  be  modeled,  as  well  the  errors  that  affect 
those  states.  However,  too  much  complexity  in  the  modeling  can  also  unnecessarily  bur¬ 
den  the  filter. 

Reference  Frames 

Since  attitude  information  is  the  goal  of  this  system,  there  must  be  some  state  representa¬ 
tion  of  vehicle  attitude.  The  vehicle  orientation  with  respect  to  the  navigation  frame  (nav 
frame)  is  of  primary  interest.  The  body  frame  will  be  represented  with  the  subscript  B 
while  the  nav  frame  will  be  denoted  with  an  N.  The  nav  frame  used  for  this  application  is 
the  Local  Vertical  Local  Horizontal  (LVLH)  frame.  This  frame  is  constructed  with  the  X 
axis  pointed  in  the  velocity  direction,  the  Z  axis  toward  nadir  and  the  Y  axis  out  of  plane. 
Two  other  reference  frames  important  in  this  study  are  the  Earth-Centered-Earth-Fixed  (E) 
frame  and  the  Inertial  (I)  reference  frame.  The  E  frame  is  the  frame  in  which  the  GPS  SV 
positions  are  reported.  It  has  its  X  axis  in  the  equatorial  plane  pointing  toward  the  Green¬ 
wich  Meridian,  its  Z  axis  aligned  approximately  with  the  Earth’s  spin  axis,  and  the  Y  axis 
completes  the  right-handed  triad  [23].  The  I  frame,  for  short  time  periods,  can  be  defined 
as  the  E  frame  minus  the  Earth’s  rotation  about  the  Z  axis.  The  following  figure  summa¬ 
rizes  the  above  descriptions.  Note  that  the  only  rotation  between  the  body  frame  and  nav 
frame  is  a  yaw  rotation. 
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Figure  3.1:  Summary  of  Reference  Frames 

Attitude  State 

As  mentioned,  the  attitude  state  is  of  fundamental  importance  in  this  study.  There  are  vari¬ 
ous  methods  of  representing  the  vehicle  body  with  respect  to  the  nav  frame.  The  most 
common  representation  is  through  the  actual  physical  angles,  known  as  Euler  angles, 
between  the  nav  and  body  frames.  The  disadvantages  of  this  method,  however,  include 
singularities  and  burdensome  trigonometric  terms.  Another  very  common  method  used  in 
many  applications  is  the  direction  cosine  matrix.  Although  it  eliminates  these  problems,  it 
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contains  nine  parameters,  some  of  which  are  redundant.  For  this  application  the  quater¬ 
nion  representation  of  the  attitude,  which  overcomes  these  disadvantages,  will  be  adopted. 

The  benefits  of  the  quaternion  attitude  representation  have  been  known  and  applied  since 
Hamilton  first  developed  it  more  than  150  years  ago.  However,  its  use  in  Kalman  filtering 
has  had  some  difficulty  due  to  the  constraint  that  the  four  quaternion  parameters  have  unit 
norm  which  results  in  the  singularity  of  the  attitude  portion  of  the  covariance  matrix.  Lef- 
ferts,  Markley  and  Shuster  were  able  to  overcome  this  problem,  and  the  use  of  quaternion 
attitude  representation  in  filtering  applications  has  become  more  commonplace  in  the  last 
15  years.  For  more  information  on  this  subject  see  [15]. 


The  quaternion,  q,  is  defined  by  four  parameters  which  represent  some  rotation  d>  about 
an  axis  of  rotation,  e : 


9  = 


, 

e  •  sin  — 


cos  — 


<ll 

<h 

^3 


<1 


(3.28) 


where  q  is  the  vector  part  of  the  quaternion,  and  q4  is  the  scalar  portion.  The  four  quater¬ 
nion  parameters  satisfy  the  aforementioned  constraint: 


2  2  2  2 
<h  +<h  +  <l3  +  ‘l4  =  1 


(3.29) 


As  stated  previously  in  the  beginning  of  the  chapter  in  section  3.2,  we  wish  the  filter  to 
keep  track  of  the  best  estimate  of  the  difference  between  the  true  attitude  and  the  estimated 
attitude.  Equivalently,  we  can  say  we  want  to  keep  track  of  the  estimation  error  between 
the  nav  and  body  frame,  i.e., 


CBN  ~  CB&CBN 


(3.30) 


where  Cfi^  can  also  be  viewed  as  the  error  in  rotation  from  nav  to  body,  or  CBn  ■  This 
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error  can  be  expressed  through  small  Euler  angle  rotations  from  the  body  as 

Cbn  ~  [1  +  [50]x]  (3.31) 

where  [80]x  is  the  vector  of  the  small  Euler  angles  in  cross-product  form: 

0  \|/  —0 

[60]x=  -i|/  o  cp  (3.32) 

9  -q>  0 

and  \)/,  0,  and  <|>  are  the  small  Euler  yaw,  pitch,  and  roll  angles,  respectively.  The  small 
angle  error  vector  is  related  to  the  quaternion  by 

5?  =  ^50  (3.33) 

where 


bq=  (q®q  ’)  (3.34) 

and  q  represents  the  true  attitude  state,  while  q  represents  the  estimated  attitude  state.  ® 
is  the  quaternion  operator  described  by 


(3.36) 


Finally,  the  attitude  state,  §qBN,  is  the  error  quaternion  from  nav  to  body  frame  expressed 
in  body  coordinates.  The  omission  of  the  scalar  part  of  the  quaternion  is  in  accordance 


with  [15]. 
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Model  Replacement  State 

As  mentioned  previously  in  section  2.3,  due  to  the  quality  of  the  low  cost  gyros,  the  tradi¬ 
tional  “model  replacement”  mode  is  not  used  in  this  application.  Instead  they  are  used  as 
any  other  measurement  with  appropriate  intensities  and  sensitivities.  This  “gyro  replace¬ 
ment”  mode  requires  the  incorporation  of  a  state  to  capture  the  modeled  dynamics.  For 
this  reason  the  dynamic  modeling  state,  the  attitude  rate  error  state,  is  included  in  the  filter 
design. 

The  angular  rate  error  is  defined  as: 

^—bn~  Qbn  ~  ®bn  (3.37) 

which  is  the  difference  between  the  true  and  estimated  rate.  This  state  will  play  an  impor¬ 
tant  role  in  state  propagation  which  will  be  needed  to  form  the  predicted  state.  The  dynam¬ 
ics  of  this  state,  largely  influenced  by  the  gravity  gradient  torque  for  space  applications, 
will  be  provided  in  section  3.3. 


Other  States 

The  rest  of  the  states  serve  to  model  the  errors  of  the  previous  chapter.  The  total  states  of 
the  filter  depend  on  the  type  of  environment  the  filter  is  used.  The  baseline  error  state  is 
omitted  for  the  roof-top  filter.  Also,  no  disturbance  torque  state  is  included  following  the 
sub-optimal  analysis  performed  by  Puri  [20],  The  only  disturbance  torque  modeled  is  the 
inclusion  of  gravity  gradient  affects  on  the  dynamics  in  the  space  filter.  The  total  state 
matrix  then  is: 


x  = 


WBN  SsC  K  A|3  AE,  (54) 


IT 


(3.38) 
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3.3  State  Dynamics 

The  dynamics  of  the  modeled  system  will  be  used  to  propagate  the  state.  First  the  non-lin¬ 
ear  dynamics  will  be  presented.  The  second  sub-section  describes  the  linearized  dynamics 
of  the  Extended  Kalman  filter. 


3.3.1  Non-Linear  Dynamics 
Attitude 


The  rate  of  change  of  a  quaternion  can  be  found  to  be: 


=  ^(fi>(0®«(0) 


(3.39) 


This  equation  can  be  used  to  express  the  attitude  kinematics  by  separating  the  angular 
velocity  into  two  components:  the  angular  velocity  of  the  body  with  respect  to  the  inertial 
frame  and  the  rate  of  the  body  with  respect  to  the  nav  frame: 


d  B  1  b  „  B  1  fi  b 

1Abn~  2~BI®  ~BN~  2-BN®  ~Nl 

which  can  be  more  conveniently  written  as  the  extension  of  equation  3.37: 

d  B  1  B  _  B 
dt-BN  ~  2 -BN®  4  BN 

due  to  the  relation 


(3.40) 


(3.41) 


—BI=  —BN  +  -NI  (3.42) 

The  dynamics  of  the  attitude  error  are  obtained  by  manipulation  of  equations  3.35,  3.38, 
and  3.40  to  yield  the  following  for  the  vector  part  of  the  quaternion  [4,  15,  25]: 

dt^BN  =  [-bn^x^Qbn*  2^-bn  (3.43) 

Attitude  Rate 

The  derivation  of  the  attitude  error  rate  dynamics  is  somewhat  more  involved.  In  the  inter- 
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est  of  clarity,  the  non-linear  equation  for  angular  rate  error  will  be  presented  without  deri¬ 
vation.  For  a  derivation  of  the  following  equations  see  Appendix  A.  The  equation  is: 

=  /‘1[6ffl2m([rx]/-[/r]x)[fx]]8q®K  (3.44) 

+  r'({Im,u]x  -  x  /])(°<4v  +  Sen®,)  -  8(d“,  x  So)“„ 

This  equation  contains  the  dynamic  modeling  that  will  aid  the  filter  in  attitude  extrapola¬ 
tion. 


3.3.2  Linearized  Dynamics 

The  non-linear  attitude  error  states 

S<  =  h™  <3-45> 

have  the  dynamics  presented  in  equations  3.43  and  3.44.  The  linearized  dynamics  matrix 
in  the  Extended  Kalman  filter  will  be  formed  by  partial  differentiation  of  the  nonlinear 
equations  with  respect  to  each  state: 


(3.46) 

(3.47) 

(3.48) 

(3.49) 

(3.50) 
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The  other  states,  which  are  the  gyro  and  GPS  error  states, 


*2=  (3.51) 

are  governed  by  their  respective  dynamics: 


The  process  noise  follows  the  same  formulation: 

wi  =  [°lx3  (w'Scd)]  (3*33) 


with 


(3.54) 


13x3  0  0  0  0 

0  i3x3  0  0  0 

g2  =  0  0  l18xl8  0  0 

®  ®  ®  l3x3  ® 

0  0  0  0  i3x3 


(3.55) 


(3.56) 
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The  reason  the  attitude  rate  error  process  noise  is  in  parenthesis  is  because  the  noise  is 
only  added  for  the  space  filter.  The  parenthesis  around  the  baseline  length  error  process 
noise  serves  to  indicate  that  this  state  is  also  only  modeled  in  the  space  filter.  It  is  assumed 
that  the  dynamics  of  the  roof-top  filter  are  fairly  well  known,  whereas  mismodeling  orbital 
dynamics  is  much  more  likely,  especially  in  the  presence  of  disturbance  torques.  Finally, 
the  total  dynamics  are: 


jci  =  Fjic  +  GjiV] 

(3.57) 

X2  =  F2x2  +  G2w2 

(3.58) 

3.4  Measurement  Prediction  and  Sensitivities 


The  key  to  maximizing  the  benefit  that  measurements  provide  to  attitude  estimation  in 
Kalman  filtering  lies  in  the  fidelity  of  the  measurement  model.  Better  measurement  mod¬ 
els  will  provide  better  sensitivities  which  will  result  in  more  optimal  filter  gains  to  bring 
the  state  closer  to  truth.  Recall  that  the  optimal  gain  correction  applied  by  the  filter  is 
based  on  the  measurement  sensitivity  matrix: 

dh\ 


(3.14) 


The  first  sub-section  will  present  the  GPS  measurement  prediction  and  sensitivities.  The 
second  section  will  derive  the  gyro  measurement  prediction  and  sensitivities. 


3.4.1  GPS  Measurement  Prediction  and  Sensitivities 

The  GPS  measurement  prediction  and  sensitivities  have  been  described  in  detail  in  [20], 
They  are  summarized  here  for  convenience.  To  create  the  GPS  measurement  residual, 

h0PS  =  A9(f*)-A9(f*_  )  (3.59) 

the  predicted  phase  is  formulated  by  applying  the  following  equation: 
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A<fc,;  =  (<)T[pLJ-A%  +  «v 


(3.60) 


The  integer  ambiguity,  Akjj ,  is  predicted  with  the  integer  ambiguity  resolver  mentioned  in 
Chapter  1.  The  predicted  error,  e^,  and  the  estimated  LOS  vector,  est ,  are  formulated 
in  [20]  but  their  definitions  are  summarized  with  the  following  two  equations: 


a  =  xa^.«].0+Afc  +  Aty+(«L[P®.„])8*l 


(3.61) 


r i  B~N 

tPy,  est^  ~  G  ’  P  j  (3.62) 

X  is  the  calibration  from  the  bias  map  described  in  Chapter  2.  Following  the  derivation  in 
[20],  the  sensitivity  of  the  GPS  measurement  residual  to  the  attitude  state  is: 


dz;  : 


-.jfV  =  2-Wj,esth-bi0 
d8qBN 


(3.63) 


There  is  no  phase  measurement  sensitivity  to  the  angular  rate  or  gyro  bias.  The  sensitivi¬ 
ties  for  the  GPS  error  states  were  given  in  Chapter  2  and  will  not  be  repeated  here. 


3.4.2  Gyro  Measurement  Prediction  and  Sensitivities 

The  following  section  describes  how  the  low  cost  gyro  measurements  are  treated,  and  the 
prediction  and  sensitivities  are  derived. 


Prediction  and  Sensitivities 

The  gyro  measurement  residual  is  simply  the  actual  measurement  minus  the  predicted 
measurement,  or: 

~zg  =  &M-QM  (3.64) 

where 


~  @BI  +  bg+X>g 


(3.65) 
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(3.66) 


so  by  applying  relation  3.42: 

zg  =  (QBbn  -  &Bbn)  +  (UBNJ  -  foBN!)  +  bg  +  vg  (3.67) 

B 

The  inertial  to  nav  predicted  rate,©w ,  can  be  found  by  realizing  that 

4  =  Cbn(qnni)  =  [1  +  [80]x](  to")  (3.68) 

Application  of  the  above  relations  along  with  equation  3.32,  yields: 


zg  -  +  [26^rBAr]x  •  ((QNl)  +  bg  +  t)g 

This  can  be  differentiated  to  find  the  measurement  sensitivities: 

(3.69) 

Js  =-2[S&lx 

36qBA, 

(3.70) 

dz„ 

B  “  (!)3x3 

O0(0BN 

(3.71) 

di‘  -  m 

3is  (1)j" 

(3.72) 

3.5  Filter  Algorithms 

The  Extended  Kalman  filter  described  in  the  previous  sections  is  used  in  two  scenarios. 
The  first,  the  roof-top  filter,  will  test  the  integration  issues  on  a  static  mount  on  the  roof  of 
Draper  Laboratory.  The  second  scenario  will  study  how  these  issues  are  affected  by  the 
high  dynamic  environment  of  a  LEO  by  integrating  the  filter  with  a  space  simulation.  The 
integration  issues,  such  as  power  consumption  and  real  MM  Gyro  data  integration,  are  the 
focus  of  this  study  and  will  be  tested  on  the  two  types  of  filters.  The  basic  filter  algorithm 
is  presented  in  the  following  figure. 


56 


Initialization: 


Figure  3.2:  Basic  Filter  Algorithm 
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3.5.1  Roof-top  Filter 

The  roof-top  filter  is  used  to  post-process  GPS  data  taken  from  the  Trimble  TANS  Vector 
receiver  on  the  roof  of  Draper  Laboratory.  Post-processing  is  required  for  a  number  of  rea¬ 
sons.  First,  the  MM  Gyros  are  still  not  ready  for  real-time  integration  with  other  instru¬ 
ments  and  the  equipment  required  to  run  the  gyros  could  not  be  taken  up  to  the  roof.  Also, 
the  temperature  fluctuations  on  the  roof  would  add  an  additional  degree  of  uncertainty  in 
the  gyro  measurements  that  would  preclude  any  useful  analysis.  Post-processing  is  also  a 
necessity  because  of  the  receiver  used.  The  Trimble  is  designed  to  calculate  a  point  solu¬ 
tion  from  each  differential  phase  measurement.  This  means  that  it  has  its  own  algorithms 
that  have  nothing  to  do  with  the  Kalman  filter  described  above.  Fortunately,  it  does  output 
the  raw  differential  phase  measurement  which  the  IGPS/MM  Gyro  ADS  uses  as  one  of  its 
two  measurement  sources.  Since  the  receiver  outputs  these  measurements  at  1  Hz,  post¬ 
processing  is  required  to  vary  the  frequency  of  GPS  measurement  use. 

The  roof-top  filter  lacks  the  baseline  estimation  state.  Temperature  fluctuations  on  the  roof 
are  small  enough  not  to  affect  the  size  of  the  rigid  aluminum  array  on  which  the  antennas 
are  mounted.  Since  there  is  more  multipath  on  the  roof  (the  shed,  heat  stacks,  and  sky  light 
mentioned  in  Chapter  2),  analysis  of  data  takes  from  the  roof  have  been  done  with  and 
without  an  extra  multipath  state.  This  will  be  described  later. 

3.5.2  Space  Simulation 

In  order  to  validate  this  system  for  use  in  space,  Puri  developed  an  orbital  simulation  that 
incorporates  the  filter  in  a  LEO  environment  [20].  This  same  simulation  (with  minor 
changes)  will  be  used  to  test  the  integration  issues  which  are  the  focus  of  this  study. 
Importantly,  it  will  also  verify  the  ability  of  the  IGPS/MM  Gyro  ADS  to  handle  high  fre¬ 
quency  dynamics:  something  the  roof-top  version  cannot  do.  For  more  information  on  the 
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details  and  mechanics  of  the  orbital  simulation,  see  reference  [20].  A  summary  of  the  sim¬ 
ulation  is  presented  below  for  convenience. 


Simulation  Parameters 


Performance  Statistics 


Figure  3.3:  Summary  of  Orbital  Simulation 

GPS  measurement  errors  for  the  simulation  are  generated  using  the  dynamic  error  models 
developed  in  Chapter  2.  The  orbit  used  in  this  study  will  be  a  typical  LEO  used  by  small 
communications  satellites.  The  orbit  is  assumed  to  be  circular  with  variable  inclination, 
semi-major  axis,  and  longitude  of  ascending  node.  No  orbit  pertubations  are  modeled.  For 
the  space  simulations  the  inclination  chosen  is  89  degrees,  and  the  altitude  is  421  nm, 
which  corresponds  to  the  orbit  of  the  Iridium  satellite  constellation.  The  GPS  constellation 
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is  propagated  using  the  ICD-GPS-200  almanac  equations  [20], 


60 


Chapter  4 


Covariance  Analysis 

4.1  Overview 

This  chapter  presents  a  covariance  analysis  of  the  Kalman  filter  to  analyze  different  inte¬ 
gration  schemes  for  the  ADS.  One  area  of  concern  is  the  effect  of  using  varying  GPS  mea¬ 
surement  update  rates.  Because  the  goal  of  this  application  is  to  provide  an  ADS  of 
reasonable  quality  and  low  power,  the  trade-offs  between  performance  and  power  usage 
must  be  analyzed.  Power  consumption  will  be  addressed  Chapter  7,  but  a  baseline  assump¬ 
tion  of  the  covariance  analysis  is  that  most  of  the  power  used  by  the  IGPS/MM  Gyro  ADS 
is  consumed  by  the  interferometric  GPS  receiver.  This  assumption  seems  valid  knowing 
that  the  TANS  Vector  receiver  itself  uses  four  Watts  at  its  12-18  Volt  operating  range, 
while  the  gyro  system  will  use  well  under  one  Watt  of  power.  Covariance  analysis  will  be 
done  on  the  roof-top  filter  as  well  as  on  the  space  filter.  This  analysis  will  begin  by  looking 
at  the  covariance  matrix  of  the  EKF  when  run  with  real  measurements;  however,  a  linear¬ 
ized  covariance  analysis  will  be  conducted  for  most  of  the  study.  The  linearized  analysis  is 
performed  because  of  the  great  computational  burden  required  for  Monte  Carlo  runs  with 
the  EKF. 

4.2  Methodology 

The  following  sub-sections  describe  the  manner  in  which  the  covariance  analysis  in  this 
study  was  conducted. 
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4.2.1  Equations 

The  Kalman  filter  can  be  used  not  only  as  an  estimator,  but  also  as  a  means  to  conduct  sys¬ 
tem  error  analyses  to  compare  different  plausible  designs.  The  advantage  of  a  Kalman  fil¬ 
ter  covariance  analysis  is  that  the  estimate,  x ,  does  not  have  to  be  propagated  along  with 
the  covariance  matrix,  P.  To  obtain  the  steady  state  value  of  the  covariance  we  simply 
apply  the  recursive  equations  of  the  previous  chapter  [1].  The  following  figure  demon¬ 
strates  the  recursive  loop  for  propagating  the  error  covariance  matrix.  Note  that  steady 
state  analysis  techniques,  such  as  solving  the  closed  loop  Lyapunov  equation,  will  not 
work  in  this  case.  This  recursive  loop  must  be  used  due  to  the  different  frequencies  at 
which  the  two  discrete  measurements  are  input  to  the  filter. 


Figure  4.1:  Recursive  Covariance  Propagation 
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4.2.2  Linearization 


In  the  EKF  the  non-linear  system  is  linearized  about  a  trajectory  that  is  continuously 
updated  with  each  state  estimate  obtained  from  the  measurements.  In  the  above  algorithm, 
however,  the  state  estimate  is  not  propagated  and  there  are  no  measurements.  To  linearize 
the  system,  a  nominal  point  in  the  trajectory  is  chosen  that  is  representative  of  what  the  fil¬ 
ter  would  experience  if  it  was  running  in  real-time  using  real  measurements.  Since  the 
non-linearity  of  the  system  is  caused  by  the  GPS  SV’s  motion  through  the  visible  sky,  a 
frozen  nominal  geometry  is  chosen.  It  consists  of  six  SV’s  with  an  Attitude  Dilution  of 
Precision  (ADOP)  of  about  .053  (see  Appendix  C).  A  factor  of  safety  is  included  by  inval¬ 
idating  the  lowest  satellite  in  elevation.  This  nominal  geometry  is  shown  in  the  sky  plot 
below: 


180 


Figure  4.2:  Nominal  Geometry  Used  in  Linearization 

The  sky  plot  is  a  representation  of  GPS  SV’s  as  seen  by  the  receiver.  Elevation  is  measured 
from  0°  to  90°  where  0°  is  the  horizon.  Azimuth  is  measured  in  standard  north,  south,  east 
and  west  where  360°  is  north.  All  the  measurement  sensitivities  for  the  linearized  covari- 
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ance  analysis  will  be  derived  from  this  geometry. 


Verifying  the  Linearization 

To  ensure  that  the  above  linearization  technique  will  yield  accurate  variances,  it  must  be 
compared  to  actual  EKF  variances.  The  following  plot  shows  the  results  of  the  roof-top  fil¬ 
ter  for  a  standard  run.  The  differential  phase  data  was  recorded  on  the  Draper  Lab  roof  on 
January  14,  1998,  at  a  rate  of  1  Hz  from  the  Tans  Vector.  The  IGPS  antenna  array  was  the 
0.6  meter  square  mentioned  in  the  first  chapter.  Simulated  gyro  measurements  were  inte¬ 
grated  every  2  Hz.  The  gyro  measurements  were  produced  according  to  the  gyro  measure¬ 
ment  model  presented  in  Chapter  2.  Unless  otherwise  stated,  the  following  tests  are  done 
assuming  a  MM  Gyro  bias  of  5°/hr  and  an  ARW  of  .25°/hr1/2.  The  EKF  la  standard  devi¬ 
ation  lines  are  overlaid  about  the  mean  of  the  solution.  The  mean  had  to  be  used  due  to  the 
uncertainty  of  the  actual  truth  orientation  of  the  array,  and  it  is  computed  by  taking  the 
average  value  of  the  solution  over  the  entire  interval.  The  run  comprises  approximately 
three  and  a  half  hours  of  data.  The  satellite  geometries,  measurement  validities  and 
ADOPs  for  this  run  are  given  in  Appendix  B. 
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KF  Attitude  Solution  for  Entire  Data  Take 


Figure  4.3:  EKF  Results  for  01/14/98  GPS  Data 

In  order  to  see  the  effect  of  reducing  the  GPS  measurement  update  frequency  on  the  vari¬ 
ance  of  the  solution,  various  runs  on  the  same  data  were  performed  with  the  only  differ¬ 
ence  being  the  frequency  of  GPS  measurement  use.  In  other  words,  instead  of  using  all  the 
differential  phase  data  recorded,  only  an  amount  commensurate  with  the  GPS  measure¬ 
ment  update  frequency  was  used.  This  was  possible  since  the  differential  phase  data  is 
time-tagged  by  the  receiver.  The  next  figure  summarizes  the  effect  this  had  on  the  EKF 
attitude  error.  The  steady  state  value  was  obtained  by  averaging  the  value  of  the  attitude 
variance  during  the  last  few  seconds  of  the  run. 
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EICF  Variances  for  7  Different  GPS  Frequency  Runs  (Hz) 


Figure  4.4:  EKF  Steady  State  Standard  Deviations  for  Varying  GPS  Use 

We  see  that  attitude  errors  did  not  increase  much  from  using  GPS  measurements  once 
every  second  to  once  every  minute.  Using  even  less  GPS  updates  than  the  1/64  Hz  rate, 
however,  does  cause  a  filter  reset  (see  Figure  3.2).  We  would  expect  the  linearized  covari¬ 
ance  analysis  to  show  similar  results.  The  following  figure  displays  the  steady  state  results 
of  the  linearized  recursive  covariance  analysis: 
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Steady  State  Variances  for  7  Different  GPS  Frequency  Runs  (Hz) 


Figure  4.5:  Steady  State  Error  of  Linear  Analysis 

From  this  figure  we  conclude  that  the  linearized  covariance  results  are  indeed  credible 
since  the  general  trend  and  magnitude  of  the  standard  deviations  are  similar.  Modest  dis¬ 
crepancies,  as  those  seen  in  the  yaw  values,  are  the  result  of  using  the  same  geometry  ver¬ 
sus  a  changing  environment.  The  yaw  performance  is  always  better  than  roll  or  pitch  since 
the  IGPS  array  is  orthogonal  to  the  yaw  axis. 

We  would  expect  that  continually  reducing  the  frequency  of  GPS  measurements  would 
indeed  degrade  the  accuracy  of  the  attitude  solution.  This  is  confirmed  by  the  following 
plot  which  carries  the  trend  of  reducing  GPS  use  still  further: 
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Steady  State  Variances  with  Decreasing  GPS  Frequency  (Hz) 


Figure  4.6:  Linear  Analysis  for  Continually  Less  GPS 

The  results  of  this  graph  serve  to  show  that  for  the  region  in  which  the  EKF  will  not  reset 
(any  GPS  frequency  use  below  1/64  Hz),  the  GPS  update  frequency  has  negligible  effect 
on  the  error  bounds.  A  filter  reset  is  not  necessarily  a  bad  thing,  but  it  does  prohibit  steady 
state  variance  analysis  since  the  covariance  matrix  is  re-initialized.  Also,  the  aim  of  this 
analysis  is  to  evaluate  the  impact  of  mitigating  the  use  of  the  “power  hungry”  GPS  inter¬ 
ferometer  on  performance.  Resetting  the  filter  calls  the  integer  ambiguity  resolver  which 
requires  an  additional,  although  slight,  computational  burden. 


The  above  results  were  for  the  static  0.6  meter  per  side  square  platform.  Similar  analysis 
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can  be  done  with  the  space  filter  in  order  to  determine  the  effects  of  dynamics  on  the  error. 
For  this  analysis,  the  space  filter  will  be  placed  in  the  LEO  orbit  described  in  Chapter  3. 
The  results  of  one  run  of  7000  seconds  are  shown.  This  particular  simulation  run  used  the 
0.6  meter  IGPS  array. 
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Figure  4.7:  EKF  Orbital  Simulation  Results 

For  the  above  run  the  GPS  update  frequency  was  1Hz.  The  1  a  standard  deviations  over¬ 
lays  are  included  along  with  the  errors.  The  errors  are  derived  by  subtracting  the  estimate 


Std.  Dev.  in  Degrees 


from  the  truth.  The  process  noise  intensity  used  on  the  attitude  rate  was  1.2  X  10'5  N-m. 
When  GPS  measurements  are  varied  from  1Hz  to  l/80Hz,  the  errors  are  bounded  in  a  sim¬ 
ilar  fashion  as  in  the  static  roof-top  case: 


KF  Steady  State  Standard  Deviation  for  Varying  GPS  Update  Frequencies 


1  1/2  1/5  1/10  1/20  1/40  1/80 

EKF  Variances  for  7  Different  GPS  Frequency  Runs  (Hz) 


Figure  4.8:  EKF  Standard  Deviations  for  Varying  GPS  Use  Sim  Runs 

The  linearized  version  of  the  above  analysis  conforms  to  the  expected  trends  as  shown  in 
the  following  plot: 
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Linear  Covariance  Analysis  for  Varying  GPS  Update  Frequencies 


Figure  4.9:  Linearized  Analysis  for  Space  Sim 

The  previous  plots  support  the  accuracy  of  the  linearized  covariance  analysis  method.  This 
linear  analysis  is  the  basis  for  the  rest  of  the  covariance  analysis  which  follows. 


4.3  Analysis 

The  conclusion  that  can  be  drawn  from  the  space  simulation  results  is  important  in  under¬ 
standing  the  performance  limits  of  the  system.  As  compared  to  the  roof  tests,  we  would 
expect  the  errors  of  the  high  dynamic  space  simulation  environment  to  increase  due  to  the 
higher  uncertainty  of  orbital  dynamics  represented  by  the  added  process  noise.  The  orbital 
errors  were  actually  a  bit  smaller  since  the  error  bounds  had  to  be  increased  on  the  roof  fil¬ 
ter  to  make  up  for  multipath  effects  on  the  roof.  What  is  important  to  note,  however,  is  that 
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the  space  filter’s  errors  were  very  comparable  to  the  roof  filter’s  errors. 

4.3.1  Unmodeled  Dynamics 

There  are  two  reasons  why  the  space  filter  errors  do  not  differ  greatly  from  the  roof-top  fil¬ 
ter  errors.  The  body-frame  motion  experienced  by  the  static  roof-top  filter  due  to  the  mov¬ 
ing  satellite  geometries  is  comparable  to  the  slow  body-frame  dynamics  of  a  spacecraft. 
As  long  as  LOS  vectors  are  updated  quickly  enough  in  the  space  simulation,  the  effects  of 
the  gravity  gradient  torque  are  not  that  large.  Another  possible  explanation  is  that  the  EKF 
in  the  space  simulation  is  based  on  a  model  of  the  dynamics  that  is  rather  simple.  If  the 
EKF  is  depending  on  the  model  too  much  to  propagate  the  state,  in  other  words  the  pro¬ 
cess  noise  on  angular  rate  error  is  low,  it  is  difficult  to  see  just  how  much  of  an  impact 
varying  GPS  measurement  frequency  or  gyro  quality  has  on  the  performance  of  the  ADS. 

The  nominal  design  process  noise  on  the  angular  rate  error,  a5(0,  of  1.2  X  10'5  N-m  inten¬ 
sity  is  an  optimistic  estimate  of  likely  orbital  conditions.  If  typical  sources  of  unmodeled 
torques  are  summed,  such  as  aerodynamic  drag,  solar  radiation  pressure,  magnetic  field 
torque,  and  internal  spacecraft  dynamics,  the  intensity  of  the  unmodeled  disturbance 

torques  might  reach  as  high  as  10'3  N-m  [26].  This  value  does  not  include  any  perturba¬ 
tions  of  the  “known”  spacecraft  attitude  control  system  torques.  In  order  to  observe  the 
effects  of  varying  the  design  process  noise,  analysis  is  done  on  a  0.6  meter  and  1  m  IGPS 

array  by  using  a  suite  of  lower  quality  gyros  (b  =  100°/hr,  x  =  2  hrs  and  ARW  =  .75°/hr17 

5  6 

2)  and  a  1  Hz  GPS  update  frequency  rate: 
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Figure  4.10:  Attitude  and  Rate  Error  as  a  Function  of  Design  Process  Noise  Intensity 

for  0.6  m  Array 


Figure  4.11:  Attitude  and  Rate  Error  as  a  Function  of  Design  Process  Noise  Intensity 

for  1  m  Array 

Attitude  error  drops  as  expected  when  the  dynamic  model  is  more  precise.  We  see,  how- 
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ever,  that  even  though  the  angular  rate  error  increases  with  increasing  angular  rate  process 
noise,  there  is  an  upper  bound  on  the  RSS  attitude  error.  This  bound  is  provided  by  the 
GPS  phase  measurements.  As  will  be  shown  in  the  next  few  pages,  having  a  better  gyro 
will  help  tighten  these  error  bounds  since  the  unknown  disturbances  will  be  better  esti¬ 
mated. 

The  error,  as  a  function  of  decreasing  GPS  updates  and  increasing  design  process  noise, 
will  now  be  examined.  The  following  3-D  bar  graph  summarizes  these  effects  for  a  0.6 
meter  IGPS  antenna  array. 


Figure  4.12:  S/C  Attitude  Error  as  a  Function  of  GPS  Update  Interval  and  Dynamic 
Model  Quality/Design  Process  Noise  for  a  0.6  m  IGPS  Array 

The  bound  that  the  GPS  measurements  had  imposed  in  Figure  4.10  is  now  stretched  with 

decreasing  GPS  updates.  The  unknown  dynamics  are  allowed  to  propagate  for  a  longer 

time,  hence  decreasing  the  performance  of  the  estimator.  This  bound,  however,  can  be 

reduced  with  a  suite  of  better  gyros.  Figure  4.13  demonstrates  the  performance  improve- 
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ments  obtained  when  increasingly  better  gyros  are  applied.  The  different  gyro  suites  tested 
are  listed  in  the  following  table. 


Name 

Bias 

(deg/hr) 

Time 

Constant 

(hr) 

ARW 

deg/hr172 

Suite  1 

1 

8 

.1 

Suite  2 

5 

8 

.25 

Suite  3 

10 

8 

.25 

Suite  4 

20 

8 

.5 

Suite  5 

40 

6 

.75 

Suite  6 

80 

4 

.75 

Suite  7 

100 

2 

1 

Table  4.1:  Description  of  Various  Gyro  Suites  Analyzed 


The  attitude  rate  process  noise  was  held  at  1.5  X  10'4  N-m  which  seems  to  be  a  reasonable 
value  that  can  be  expected  from  external  spacecraft  torques  (not  to  include  ACS  torques) 
[26]. 
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Figure  4.13:  Attitude  Error  as  a  Function  of  GPS  Update  Interval  and  Gyro  Quality 

for  a  0.6  m  IGPS  Array 

The  higher  the  gyro  quality,  the  better  the  ADS  can  track  high  frequency  dynamics  during 
GPS  outages.  Performance  and  cost  are  synonymous  with  trade-offs  in  engineering,  and 
this  system  is  no  different.  What  remains  to  be  analyzed  is  how  much  power  the  IGPS 
receiver  requires  and  how  compatible  this  cost  is  with  available/obtainable  gyro  qualities. 
This  will  be  discussed  in  Chapter  7. 

A  repeat  of  the  analysis  done  in  Figure  4.12  and  4.13  is  now  shown  using  a  1  meter  IGPS 
antenna  array.  Increasing  the  size  of  the  array  is  one  method  to  improve  this  ADS’s  perfor¬ 
mance.  This  comes  at  the  expense  of  larger  size  and  a  greater  computational  burden  to 
resolve  the  integer  ambiguities.  This  also  will  be  discussed  in  Chapter  7. 
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Figure  4.14:  S/C  Attitude  Error  as  a  Function  of  GPS  Update  Interval  and  Dynamic  Model 
Quality/Design  Process  Noise  for  a  1  m  IGPS  Array 


Figure  4.15:  Attitude  Error  as  a  Function  of  GPS  Update  Interval  and  Gyro  Quality  for  a  1 

m  IGPS  Array 
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Chapter  5 


Robustness  to  Failure 

One  important  aspect  of  a  system,  especially  when  comparing  it  to  other  viable  alterna¬ 
tives,  is  its  robustness  to  failure.  The  performance  of  a  system  must  be  analyzed  vis-a-vis 
the  redundancy  of  the  system.  This  chapter  will  look  at  some  inherent  redundancies  in  the 
IGPS/MM  Gyro  ADS  which  make  it  robust  to  failure. 

5.1  Introduction 

Most  aerospace  systems  are  designed  with  some  type  of  robustness  to  failure.  The  failure 
of  only  one  system,  such  as  the  ADS,  could  result  in  loss  of  mission  capabilities.  To  pro¬ 
vide  robustness  to  failure  most  aerospace  ADS’s  come  with  either  a  variety  of  instruments, 
such  as  a  suite  of  sun/horizon  sensors,  gyros,  magnetometers  and  star  trackers;  or  a  redun¬ 
dancy  in  the  number  of  certain  sensors,  such  as  multiple  gyros.  If  an  orthogonal  gyro  triad 
is  measuring  rate  about  three  axes,  then  a  loss  of  one  gyro  would  be  a  single  point  of  fail¬ 
ure.  Solutions  to  this  problem  can  either  be  to  provide  back-up  gyros  or  to  align  five  or  six 
gyros  in  such  a  manner  so  as  each  senses  components  of  rate  from  multiple  axes.  In  either 
case,  the  cost  of  providing  redundancy  is  higher  than  having  a  bare-bones  system. 

5.2  Advantages  of  the  T-C  Model  Replacement  Mode 

The  tightly-coupled  model  based  Kalman  filtering  used  for  this  application  actually  gives 
the  ADS  a  certain  degree  of  robustness  to  failure  not  inherent  with  traditional  model 
replacement  mode  systems.  The  following  figure  is  a  visualization  of  the  way  measure¬ 
ments  are  treated  in  this  application. 
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Figure  5.1:  Tightly  Coupled  vs.  Loosely  Coupled  ADS 

In  the  loosely  coupled  system  the  filter  is  used  to  optimally  blend  two  attitude  estimates 
previously  computed  by  the  gyros  and  IGPS.  The  IGPS  solution  is  then  a  least  squares 
point  solution.  One  disadvantage  with  this  method,  as  well  as  with  the  traditional  model 
replacement  mode,  is  that  each  strapdown  gyro  is  a  single  point  of  failure.  If  these  systems 
are  to  be  made  robust,  then  redundant  gyros  have  to  be  used  resulting  in  an  increase  in  cost 
of  the  entire  ADS.  The  tightly  coupled  integration  proposed  for  this  application  circum¬ 
vents  this  disadvantage.  Each  strapdown  gyro  output  is  treated  as  basically  one  measure¬ 
ment  out  of  the  total  number  of  measurements,  which  for  this  application  would  be  one  out 
of  twenty-one  measurements:  three  gyros  and  eighteen  phase  measurements  (six  channels 
times  three  baselines). 

5.2.1  Robustness  to  Gyro  Failure 

Since  there  is  a  dynamic  model  to  aid  in  state  propagation,  we  would  expect  continued 
operation  (albeit  degraded)  from  this  system  in  the  event  of  a  gyro  malfunction.  This 
would  be  impossible  in  a  model  replacement  mode  system  that  depended  on  the  gyro  triad 
output  for  state  propagation.  The  amount  of  performance  degradation  in  the  event  of  a 
gyro  malfunction  is  a  function  of  the  quality  of  the  MM  Gyros  and  the  IGPS  measurement 
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update  interval.  The  gyro  triad  characterized  by  Suite  2  in  Table  4. 1  is  a  good  approxima¬ 
tion  of  the  quality  of  MM  Gyros  that  may  be  achieved  in  the  near  future.  The  following 
table  describes  the  three  operating  states  that  will  be  examined  using  the  Suite  2  triad. 


Gyro  State 

Description 

Normal 

[111] 

Pitch  Out 

[101] 

Pitch  &  Yaw  Out 

[100] 

Table  5.1:  Descriptions  of  Operating  States 


The  effect  on  performance  of  using  these  three  different  operating  scenarios  are  listed  in 
the  following  table  as  a  function  of  IGPS  update  interval.  The  results  are  for  a  0.6  m2 
IGPS  array. 


IGPS 

Tlnrlotp 

Normal  State  [111] 

Pitch  Out  [  1 0  1  ] 

P  &  Y  Out  [  1 0  0  ] 

Freq 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

Roll 

Pitch 

Yaw 

1  Hz 

.2177 

.1904 

.1028 

.2177 

.2755 

.1031 

.2179 

.2774 

.1701 

1/2  Hz 

.2187 

.1912 

.1041 

.2187 

.2802 

.1043 

.2189 

.2823 

.1793 

1/5  Hz 

.2208 

.1931 

.1066 

.2208 

.2865 

.1068 

.2209 

.2887 

.1925 

1/10  Hz 

.2246 

.1966 

.1106 

.2246 

.3008 

.1108 

.2247 

.3030 

.2161 

1/20  Hz 

.2319 

.2032 

.1176 

.2319 

.3290 

.1178 

.2320 

.3313 

.2587 

1/40  Hz 

.2459 

.2155 

.1296 

.2459 

.3862 

.1297 

.2460 

.3886 

.3424 

1/80  Hz 

.2748 

.2405 

.1515 

.2748 

.5397 

.1517 

.2749 

.5424 

.5752 

Table  5.2:  Attitude  Error  in  Deg.  of  Three  Gyro  Conditions 
as  a  Function  of  Varying  IGPS  Update  Intervals 


These  statistics  can  be  better  visualized  in  the  following  figure. 
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Gyro  Suite  Status 


Figure  5.2:  Attitude  Error  in  Deg.  of  Three  Gyro  Conditions 
as  a  Function  of  Varying  IGPS  Update  Intervals 


The  first  bar  chart  in  Figure  5.2  conforms  with  values  in  Figure  4.8.  In  the  nominal  config¬ 
uration,  the  diagonal  baseline  of  the  IGPS  array  is  aligned  with  the  roll  axis  as  demon¬ 
strated  by  the  following  figure. 
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Patch  ^  B 


Figure  5.3:  Nominal  Configuration  of  IGPS  Array 

Notice  that  it  is  not  possible  to  observe  roll  about  a  baseline  using  measurements  solely 
from  that  baseline.  Since  the  large  diagonal  baseline  is  used  as  the  body  frame  roll  axis, 
roll  performance  is  worse  than  that  of  pitch  and  yaw.  (If  better  roll  performance  is  desired, 
rotation  of  the  array  to  align  the  pitch  axis  with  the  long  baseline  is  possible.  Improved  roll 
performance,  however,  would  come  at  the  expense  of  pitch  performance.)  Recall  that  yaw 
performance  is  better  because  the  IGPS  array  is  orthogonal  to  the  yaw  axis.  In  the  second 
sub-plot  of  Figure  5.2,  having  the  pitch  gyro  out  significantly  degrades  pitch  knowledge. 
With  decreased  IGPS  updates,  the  pitch  performance  worsens  faster  than  that  of  roll  and 
yaw.  The  third  plot  in  Figure  5.2  has  the  pitch  and  yaw  gyros  off-line.  Although  there  is  a 
noticeable  degradation  in  both  pitch  and  yaw  performance,  the  orthogonality  advantage 
mentioned  above  constrains  the  yaw  errors  at  high  IGPS  updates.  At  lower  IGPS  updates, 
however,  this  advantage  is  clearly  lost.  Yaw  knowledge  degradation  actually  surpasses  that 
of  pitch  in  the  lowest  IGPS  update  frequency  case.  This  is  due  to  the  coupling  effect  of  the 
dynamics. 

Although  the  loss  of  one  or  two  gyros  resulted  in  degraded  performance,  the  ADS  was  still 
able  to  provide  attitude  information.  To  quantify  the  robustness  to  failure  that  this  system 
provides,  nine  scenarios  were  created  for  this  ADS  each  having  a  different  operating  sta- 
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tus.  The  nine  cases  are  described  in  the  table  below.  Suite  2  is  the  same  gyro  suite 
described  in  Table  4. 1 . 


Case 

Description  of  Sensor  Status 

Case  1 

Normal  Suite  2,  IGPS  Update  of  1  Hz 

Case  2 

Normal  Suite  2,  IGPS  Update  of  1/10  Hz 

Case  3 

Suite  2  Pitch  Out,  IGPS  Update  of  1  Hz 

Case  4 

Suite  2  Pitch  Out,  IGPS  Update  of  1/10  Hz 

Case  5 

Suite  2  Pitch  &  Yaw  Out,  IGPS  Update  of  1  Hz 

Case  6 

Suite  2  Pitch  &  Yaw  Out,  IGPS  Update  of  1/10  Hz 

Case  7 

No  Gyros,  IGPS  Update  of  1  Hz 

Case  8 

No  Gyros,  IGPS  Update  of  1/10  Hz 

Case  9 

Suite  2,  No  IGPS 

Table  5.3:  Description  of  Operating  Scenarios 


The  RSS  steady  state  attitude  errors  associated  with  the  above  cases  are  shown  below.  All 
values  assume  an  IGPS  array  of  1  m2. 
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Figure  5.4:  Performance  Results  of  Various  Operating  Scenarios 

Some  interesting  conclusions  can  be  drawn  from  the  above  analysis.  The  IGPS  bounds  the 
errors  of  the  ADS  at  0.281  degrees  at  the  1  Hz  GPS  update  frequency  and  at  0.316  degrees 
with  a  1/10  Hz  rate.  Any  amount  of  gyro  information,  whether  from  a  full  suite  as  in  case 
1  and  2  or  just  one  roll  sensor  as  in  case  5  and  6,  adds  some  benefit  to  the  GPS  only  sce¬ 
nario.  Another  important  point  to  note,  however,  is  that  the  system  is  rendered  useless  if 
the  interferometer  malfunctions  as  shown  in  case  9.  (Value  shown  for  case  9  is  after  one 
hour.  Error  would  be  expected  to  increase  monotonically  with  time.)  GPS  outage  is  a  sin¬ 
gle  point  of  failure,  and  the  mean  time  until  failure  of  the  IGPS  sensor  will  dictate  the  sys¬ 
tem’s  failure  expectancy.  However,  even  current  off-the-shelf  GPS  systems  have  high 
failure  tolerances. 

5.2.2  Simulation  Implementation  for  Gyro  Loss 

The  space  simulation  is  utilized  in  order  to  test  the  robustness  claims  made  with  the  linear 
steady  state  covariance  analysis.  A  1/10  Hz  GPS  update  run  similar  to  that  shown  in  Fig- 
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ure  4.7  is  conducted,  only  that  this  time  the  pitch  sensing  gyro  is  rendered  useless  half¬ 
way  into  the  orbital  simulation.  The  gyro  failure  occurs  at  a  time  of  about  58  minutes  into 
the  run.  Simulation  results  are  shown  in  Figure  5.5: 


Yaw  Pitch  Roll 


0  50  100  0  50  100  0  50  100 


0  50  100  0  50  100  0  50  100 

Time  (min)  Time  (min)  Time  (min) 


Figure  5.5:  EKF  Orbital  Simulation  with  Pitch  Gyro 
Failure  at  Midpoint  (IGPS  =  1/10  Hz) 

We  see  that  when  the  gyro  fails,  the  error  bounds  on  the  rate  immediately  widen.  The  atti- 
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tude  bounds  follow  suit,  and  although  performance  is  degraded,  complete  system  failure  does  not 
occur.  In  fact,  integer  ambiguity  lock  was  still  maintained  in  this  scenario  and  no  filter  reset  was 
required. 

5.2.3  Robustness  to  Antenna  Failure 

It  was  mentioned  above  that  an  IGPS  receiver  outage  was  a  single  point  of  failure  in  this  system, 
however,  the  interferometer  will  continue  to  provide  useful  information  in  the  event  of  the  loss  of 
one  baseline.  Having  four  antennas  with  three  baselines  inherently  provides  robustness  to  an 
antenna  failure,  as  the  interferometer  needs  a  minimum  of  three  antennas  and  two  baselines  to 
operate.  The  following  figure  shows  the  loss  in  performance  when  no  measurements  are  used 
from  the  diagonal  baseline  (illustrated  in  Figure  5.3)  due  to  a  failure  in  antenna  number  two.  With 
the  loss  of  one  baseline  and  because  of  the  disqualification  of  the  invalid  SV  of  Figure  4.2,  a  max¬ 
imum  of  10  GPS  measurements  (versus  18  possible)  were  used  for  this  linear  analysis. 


Figure  5.6:  Performance  Degradation  with  Loss  of  Baseline  2 

Gyro  Suite  2  was  used  in  the  above  analysis.  Notice  the  degradation  of  performance  primarily  in 
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the  pitch  axis.  Yaw  is  also  affected,  but  roll  is  not.  As  mentioned  previously,  roll  about  a  baseline 
cannot  be  detected,  and  since  the  second  baseline  was  aligned  with  the  roll  axis,  no  roll  knowl¬ 
edge  is  lost. 

5.2.4  Simulation  Implementation  for  Baseline  Loss 

The  following  figure  contains  results  for  an  orbital  simulation  in  which  the  entire  orbit  was  run 
with  only  baselines  one  and  three  operating.  The  GPS  update  interval  was  1/10  Hz. 
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Figure  5.7:  Simulation  Results  with  loss  of  Baseline/Antenna  2 


The  results  show  that  errors,  although  slightly  degraded  from  the  nominal  array  configuration 
shown  in  Figure  4.7,  are  still  bound  and  that  no  filter  reset  was  needed. 


Chapter  6 


Real  Gyro  Data  Integration 

All  analysis  so  far  has  been  done  either  with  simulated  gyro  measurements  (from  the  mod¬ 
els  of  Chapter  2),  as  in  the  case  of  the  EKF  runs  with  GPS  measurements;  or  with  gyro 
error  parameters,  as  in  the  case  of  the  linear  covariance  analysis.  This  chapter  presents 
results  for  the  integration  of  real  MM  Gyro  data  on  the  static  roof-top  receiver. 

6.1  Post  Processing 

Since  the  Draper  Lab  MM  Gyros  are  still  in  development,  GPS  data  and  gyro  data  were 
taken  on  separate  days  and  in  separate  environments.  This  is  not  a  problem  since  the  gyro 
data  is  obtained  from  a  drift  run,  and  the  GPS  array  is  static  on  the  roof.  Dynamic  tests  can 
also  be  done  by  integrating  gyro  data  obtained  from  a  turntable  with  the  space  simulation. 
Pre-determined  rate  trajectories  can  be  fed  into  a  programmable  turntable  and  the  recorded 
gyro  outputs  can  then  be  used  as  the  gyro  measurements  in  the  space  simulation.  Only  a 
one-axis  turntable  is  needed  since  the  yaw,  pitch  and  roll  axes  can  be  done  separately  and 
with  the  same  test  instrument.  Dynamic  testing  with  real  gyro  data,  however,  could  not  be 
completed  in  time  for  the  publication  of  this  work.  Only  the  static  results  with  the  roof-top 
filter  are  presented  in  this  chapter. 

A  static  fifteen  hour  gyro  drift  test  was  temperature  compensated  and  divided  into  three 
five  hour  blocks  to  form  measurements  for  a  three  axis  gyro  suite.  The  gyro  output  was 
collected  at  every  second  and  since  the  filter  is  set  up  to  take  gyro  measurements  at  the  2 
Hz  rate,  an  interpolation  scheme  was  formed  to  fill  in  the  gaps.  This  scheme  consisted  of 
passing  the  gyro  data  through  a  lowpass  Butterworth  filter  to  separate  the  correlated  bias 
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from  the  noise.  The  interpolation  was  done  on  the  lowpass  output  and  then  corresponding 
random  noise  was  added  to  yield  gyro  data  every  half  second.  The  figure  below  shows  the 
temperature  compensated  drift  data  collected  from  a  prototype  gyro  and  its  corresponding 
correlated  and  white  noise  components. 


Residual  Gyro  Output  after  2nd  Order  Temperature  Fit 


Residual  Noise  After  Removal  of  Low  Pass  Signal 


400  500 

Time  (minutes) 


Figure  6.1:  Data  Analysis  from  a  Prototype  MM  Gyro 

The  MM  Gyro  data  had  a  bias  of  about  1  l°/hr  and  an  ARW  of  .57  deg/hr 1/2. 


6.2  Results 

Draper  Laboratory  has  made  tremendous  progress  in  the  development  of  MM  gyros  and 
will  continue  to  improve  them  from  their  present  developmental  stage.  The  results  from 
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integrating  the  above  data  into  the  roof-top  filter  are  shown  below.  The  GPS  updates  came 
once  every  ten  seconds.  Again,  although  exact  truth  was  not  known,  the  covariance  1  c 
overlays  are  centered  about  the  mean. 


KF  Attitude  Solution  for  Entire  Data  Take 


Figure  6.2:  EKF  Output  Integrating  Real  Gyro  Data  and  Real  GPS  Phase  Measure¬ 
ments  (GPS  Update  Interval  of  1/10  Hz) 

The  results  do  not  differ  significantly  from  Figure  4.3.  This  attests  not  only  to  the  impor¬ 
tance  of  the  dynamic  model  and  the  GPS  measurements,  but  also  shows  that  gyro  quali¬ 
ties,  of  the  levels  needed  for  useful  integration  for  this  ADS,  are  being  achieved.  A  run 
using  the  same  gyro  data  with  GPS  updates  at  about  once  every  minute  is  presented  next. 
In  this  figure  only  a  30  minute  frame  of  the  three  and  a  half  hours  is  shown. 
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Figure  6.3:  EKF  Output  Integrating  Real  Gyro  Data  and  Real  GPS  Phase  Measure¬ 
ments  (GPS  Update  Interval  of  1/64  Hz) 

Figure  6.3  better  illustrates  how  each  GPS  measurement  update  affects  the  estimate.  With 
this  view  we  can  see  that,  if  the  gyro  quality  continues  to  improve,  sporadic  GPS  measure¬ 
ments  need  to  be  used  only  for  calibrating  the  gyro’s  low  frequency  bias.  In  addition  to  its 
use  as  a  calibrating  device,  the  IGPS  receiver  will  provide  redundancy  to  accomodate  fail¬ 
ure  in  the  ADS.  If  the  gyro  suite  is  ever  lost,  the  ADS  could  continue  to  operate  as 
described  in  Chapter  5. 
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Chapter  7 


Power  Analysis 

7.1  Introduction 

In  order  for  the  IGPS/MM  Gyro  ADS  to  become  viable,  it  must  offer  an  efficient  power 
budget.  Chapter  4  analyzed  the  performance  of  this  system  as  GPS  use  and  gyro  quality 
were  varied.  This  chapter  attempts  to  quantify  the  actual  energy  expenditure  of  this  pro¬ 
posed  system.  Analysis  is  performed  on  simulations  of  typical  orbital  missions  that  the 
ADS  would  fly. 

7.2  Interferometer 

Although  the  differential  phase  data  collected  in  Chapter  6  was  done  with  a  Trimble  TANS 
Vector  receiver,  this  chapter  explores  the  concept  of  a  minimum  energy  interferometric 
GPS  receiver.  Guinon  characterizes  and  quantifies  the  energy  needed  for  such  a  receiver  in 
[10],  and  a  summary  will  be  presented  below.  Design  information  on  commercial  off-the- 
shelf  interferometers  like  the  TANS  can  be  found  in  [5], 

An  IGPS  will  consume  minimum  energy  if  any  receiver  subsystems  not  in  use  are  turned 
off.  Consequently,  when  GPS  update  intervals  are  increased  (as  explored  in  Chapter  4), 
less  energy  will  be  used  by  the  interferometer.  The  energy  needed  for  a  differential  phase 
measurement  is  based  on  the  energy  consumed  by  the  following  GPS  receiver  sub-sys¬ 
tems:  the  radio  frequency  electronics  (RFE),  the  digital  signal  processor  (DSP),  the  data 
dependent  processor  (DDP)  and  the  reference  oscillator.  They  are  shown  in  the  following 
figure. 
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Figure  7.1:  IGPS  Receiver  Sub-Systems 
Radio  Frequency  Electronics 

The  RFE’s,  consisting  of  low  noise  amplifiers,  one  or  more  superheterodyne  gain  stages 
and  an  analog  to  digital  converter;  consume  most  of  the  power  in  an  interferometer.  Along 
with  the  reference  oscillator  and  frequency  synthesizer,  the  power  consumption  of  the 
RFE  (Prfe)>  usin§  available  off-the-shelf  technology,  is  0.25  Watts  per  antenna.  Guinon 
shows  that  the  power  consumption  of  a  multiplexer  offsets  any  benefits  gained  by  sharing 
RFE’s. 

One  way  to  save  on  power  would  be  to  use  two  baselines  instead  of  three,  therefore  elimi¬ 
nating  the  need  for  one  antenna’s  RFE.  Phillips  shows  that  converting  the  1  m2  array  into 
one  with  only  two  baselines  and  three  antennas  is  feasible  [18].  The  following  figure  illus¬ 
trates  this  methodology. 


Figure  7.2:  Four  versus  Three  Antenna  Configuration 

There  is  only  a  small  trade-off  in  accuracy  for  losing  one  antenna  (since  the  baselines  span 
the  same  plane);  however,  with  the  power  savings  that  it  provides,  more  GPS  updates  can 
be  used.  The  increase  in  GPS  measurement  use  actually  allows  for  slightly  better  attitude 
estimation  at  equal  power  settings  with  the  three  antenna  set-up.  However,  redundancy  is 
lost  since  one  lost  antenna  would  cause  a  loss  of  three  axis  attitude  estimation  with  the 
IGPS.  (Although  one  baseline  could  still  provide  some  attitude  information,  this  analysis 
considers  less  than  three  functioning  antennas  a  single  point  of  failure  since  neither  integer 
ambiguity  nor  a  point  solution  can  be  calculated.)  Because  the  three  antenna  design  would 
lack  the  redundancy  analyzed  in  Chapter  5,  it  is  not  considered  further. 


The  energy  used  by  the  RFE’s  is  defined  by  the  product  of  power  the  hardware  consumes 
and  the  length  of  time  that  it  is  on: 


erf  ~  ”rfe  '  *  in  (7-1) 

where  TIN  is  the  length  of  the  signal  intercept.  Signal  intercept  length  is  a  function  of  the 
method  used  for  the  LI  code  alignment,  as  well  as  the  signal  to  noise  ratio,  and  is  derived 
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in  [Guinon]. 


Digital  Signal  Processor 

The  Digital  Signal  Processor  (DSP)  is  specifically  used  for  all  arithmetically  intensive 
computations  inherent  with  Doppler  delay  analysis  and  carrier  phase  measurement  of  the 
GPS  signal.  The  energy  consumed  by  this  subsystem  is 

edsp  ~  NOP'  E0P(DSP )  (7.2) 

where  NOP  is  the  number  of  operations,  or  complex  multiply  accumulates  (MAC’s) 
needed  to  compute  the  Dopple/delay  spectrum  for  all  channels  during  code  alignment, 
plus  the  number  of  MAC’s  needed  for  carrier  phase  measurement  on  all  antennas.  Opera¬ 
tions  on  the  order  of  1010  are  required  for  code  lock,  whereas  only  about  105  MAC’S  are 
needed  for  each  phase  measurement.  E0P  is  the  energy  per  operation  used  by  the  proces¬ 
sor.  Assuming  a  processor  such  as  an  Analog  Devices  ADSP-2100  DSP  Microcomputer, 
this  energy  consumption  would  be  about  4  nJ  (4  X  10'9  Joules)  per  2X8-bit  MAC.  For  a 
complete  derivation  of  the  minimum  energy  required  for  code  lock  and  phase  measure¬ 
ment  see  [10]. 

Data  Dependent  Processor 

The  Data  Dependent  Processor  (DDP)  is  used  for  all  calculations  other  than  those  per¬ 
formed  by  the  DSP.  These  include  satellite  selection,  computation  of  satellite  position 
from  a  stored  ephemeris,  estimating  the  pseudo  ranges  and  calculating  the  user  position 
and  clock  bias.  The  DDP  will  also  be  used  to  integrate  the  gyro  triad  output,  processing 
the  Extended  Kalman  filter  algorithms  and  computing  the  integer  ambiguities  when 
needed.  Although  the  digital  signal  processor  could  be  implemented  to  perform  the  DDP 
functions,  it  is  assumed  here  that  the  DDP  will  be  the  main  processor  onboard  the  satellite. 
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This  would  facilitate  integration  with  the  attitude  control  system  (ACS)  as  well  as  make  it  easier 
to  alter  the  EKF  from  ground  control  since  it  would  be  inherently  linked  with  the  Telemetry, 
Tracking  and  Command  (TT&C)  sub-system.  Another  benefit  of  having  the  DDP  be  part  of  the 
main  spacecraft  processor  is  better  line  bias  estimation  since  the  thermal  sub-system  on  the  space¬ 
craft  could  aid  in  temperature  monitoring  of  the  IGPS  cables.  The  energy  consumed  by  the  Data 
Dependent  Processor  for  running  the  IGPS/MM  Gyro  ADS  can  be  characterized  as 

&DDP  ~  N0p -  E0P(DDP )  (7.3) 

The  number  of  operations  performed  by  the  DDP  for  this  ADS  will  vary  depending  on  the  level  of 
integration  and  algorithm  complexity,  but  typical  computations  will  include  estimation  of  the 
pseudo  ranges  for  positioning  (which  is  important  for  LOS  calculation)  and  integer  ambiguity  res¬ 
olutions. 

With  past  interferometric  applications  the  integer  ambiguity  checks  imposed  a  major  computa¬ 
tional  burden.  There  are  two  reasons  why  this  will  not  be  a  problem  in  this  system.  First  of  all,  the 

DCAR  method  developed  by  Puri  only  requires  operations  on  the  order  of  104  to  resolve  the  inte¬ 
gers  and  compute  the  least  squares  point  solution  attitude  estimate  for  a  one  meter  square  array 
(which  is  the  longest  baseline  considered  here  useful  for  small  satellites).  This  is  a  great  improve¬ 
ment  over  previous  methods.  The  other  reason  has  to  do  with  the  low  probability  that  the  integer 
resolver  will  even  be  called.  The  Extended  Kalman  filter  keeps  track  of  the  ambiguities  by  form¬ 
ing  the  differential  phase  prediction.  The  only  time  the  integer  resolver  is  needed  is  when  the  mea¬ 
surement  residuals  become  large  due  to  the  phase  prediction  not  matching  the  phase 
measurement.  Situations  where  this  could  conceivably  occur  might  arise  due  to  unmodeled  con¬ 
trol  torquing  of  the  vehicle.  To  determine  the  number  of  integer  resolutions  this  system  is  likely  to 
expect,  the  aforementioned  orbital  simulation  was  used.  In  the  prior  simulations,  it  was  assumed 
that  the  control  torques  were  exactly  modeled  (at  least  to  the  level  of  the  process  noise  imposed  on 
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the  rate  error  state).  The  nominal  control  used  in  the  simulation  is  a  simple  deadband  controller 
that  keeps  the  spacecraft  aligned  within  a  certain  deadband  angle,  0DB.  Now  an  unmodeled  con¬ 
trol  is  added  within  the  orbit.  This  attitude  control  torque  is  typical  of  what  a  satellite  in  an  Irid- 
ium-type  constellation  could  expect  in  order  to  accomplish  mission  requirements  such  as  satellite- 
to-satellite  or  satellite-to-ground  node  tracking.  Figure  7.3  shows  the  true  yaw,  pitch  and  roll 
angles  and  rates  experienced  during  a  simulation  of  one  oribt. 


Yaw  FMtch 


Figure  7.3:  Simulation  Attitude  History 

Notice  the  larger  magnitudes  during  the  mission  dependent  control  maneuvers,  especially  in  pitch. 
Figure  7.4  displays  the  attitude  errors  along  with  the  1  a  standard  deviation  overlays.  The  simula¬ 
tion  was  run  at  a  2  Hz  gyro  update  frequency.  GPS  update  frequency  was  1/10  Hz  except  during 
the  mission  maneuvers  where  it  was  increased  to  1  Hz. 
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Notice  the  increase  in  the  rate  error  variance,  caused  by  an  increase  in  the  process  noise,  to 
account  for  the  unmodelled  control  torques.  Figure  7.5  and  7.6  show  the  history  of  the  rest  of  the 
state  errors. 
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Notice  the  decrease  in  the  number  of  valid  measurements  shown  in  Figure  7.5  during  the  mission 
dependent  maneuvers.  Even  with  less  available  measurements  the  EKF  was  able  to  hold  integer 
lock.  The  integer  ambiguity  resolver  was  called  twice  during  the  entire  orbit.  The  first  time  was 
the  initialization.  It  was  called  a  second  time  during  the  first  maneuver  period  when  GPS  measure¬ 
ments  were  scarce.  However,  since  the  difference  in  the  output  of  the  DCAR  was  within  a  set 
threshold  of  the  estimate,  a  filter  reset  was  not  called  and  the  variances  and  states  were  allowed  to 
keep  propagating  unchanged.  Similar  results  were  experienced  when  the  GPS  update  intervals 
were  kept  at  1/10  Hz  throughout  the  mission  dependent  maneuvers. 

As  shown  above,  the  number  of  operations  that  the  DDP  will  have  to  compute  will  not  be  signifi¬ 
cantly  increased  by  the  ambiguity  resolver.  Calculating  the  four  pseudo  ranges  (three  for  position, 
one  for  time)  and  running  the  filter  dependent  algorithms  will  be  largest  source  of  operations 
required.  To  obtain  a  safe  estimate  of  the  magnitude  of  necessary  calculations,  the  floating  point 
operations  count  from  the  above  simulation  in  Matlab  is  used.  The  simulation  required  an  average 

of  105  calculations  per  second.  It  must  be  noted  that  the  above  simulation  performed  many  more 
operations  than  would  be  required  of  this  ADS.  (These  include  operations  such  as  simulating  the 
gyro  and  phase  measurements  as  well  as  computing  environmental  and  system  truth  states).  Using 
this  number  on  a  typical  DDP,  such  as  a  Cirrus  Logic  ARM  7100FE  32-bit  RISC  processor,  can 
give  an  idea  of  the  power  consumption  required  of  this  sub-system.  The  above  device  consumes 
15  mWatts  at  an  18  MHz  processor  rate.  Allowing  32  cycles  for  one  floating  point  operation 
would  mean  less  than  3  mWatts  power  consumption.  Since  this  is  such  a  small  amount  under  such 
conservative  projections,  the  DDP  power  expenditure  required  for  this  ADS  will  be  ignored  for 
the  rest  of  this  analysis. 
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7.3  MM  Gyro 

The  MM  Gyros  currently  being  developed  at  Draper  Laboratory  are  being  packaged  together  with 
the  ASICs  that  run  the  electronics.  Current  prototypes  operate  at  roughly  200  mWatts  per  channel. 
New  ASICs  should  cut  this  power  down  by  a  factor  of  four  by  the  end  of  1998  [13]. 

7.4  Total  Power  Consumption 

From  the  above  analysis  total  power  consumption  for  the  entire  IGPS/MM  Gyro  ADS  can  be 
quantified.  The  total  power  will  be  a  combination  of  the  gyro  triad  expenditure  plus  the  power 
consumption  of  the  IGPS  at  the  desired  update  interval.  Adapting  analysis  performed  by  Guinon, 
the  following  figure  describes  the  energy  consumption  of  the  IGPS. 


Figure  7.7:  IGPS  Energy  Consumption  vs.  Time  Since  Last  Measurement 

The  following  figure  illustrates  the  total  ADS  power  consumption  at  varying  GPS  update  inter¬ 
vals.  Gyro  triad  consumption  of  150  mWatts  is  assumed. 
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From  the  covariance  analysis  performed  in  Chapter  4  (reference  Figure  4.6,  for  example),  we  see 
that  the  attitude  estimation  errors  increase  exponentially  with  decreasing  GPS  use.  Attitude  esti¬ 
mation  errors  do  not  significantly  grow,  however,  until  GPS  update  frequencies  are  lower  than  1/ 
30  Hz.  Using  a  GPS  update  interval  of  1/10  Hz  would  provide  the  performance  capabilities 
detailed  in  Chapter  4  with  a  minimal  ADS  power  consumption  of  200  mWatts. 
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Chapter  8 


Conclusions 

The  goal  of  this  thesis  has  been  to  address  some  of  the  integration  issues  associated  with 
the  IGPS/MM  Gyro  ADS.  Of  chief  concern  was  the  power  consumption  of  the  interferom¬ 
eter  and  the  robustness  of  the  system  to  component  failures.  Successful  tightly  coupled 
integration  of  real  micro-mechanical  gyro  data  with  real  GPS  phase  measurements  was 
also  performed. 

8.1  Summary  of  Results 

The  performance  capabilities  of  the  IGPS/MM  Gyro  ADS  are  a  function  of  various 
parameters.  These  include  gyro  quality,  dynamic  model  quality  and  the  amount  of  interfer¬ 
ometer  use.  The  cost  in  power  to  operate  this  system  is  a  function  of  the  amount  of  inter¬ 
ferometer  use  and  the  type  of  IGPS  receiver  (off-the-shelf  versus  minimal  energy 
receivers). 

A  linear  covariance  analysis  was  used  to  quantify  performance  capabilities  of  the  system 
with  variations  in  the  above  stated  parameters.  It  was  determined  that  for  GPS  phase  mea¬ 
surement  updates  in  the  range  of  one  every  second  to  one  every  30  seconds,  the  attitude 
performance  variation  was  less  than  .05  degrees.  This  justifies  the  development  of  a  mini¬ 
mum  energy  receiver  to  take  advantage  of  the  potential  power  savings  of  less  frequent 
GPS  measurements. 

Effects  of  varying  the  model  quality  of  the  EKF  were  also  quantified.  It  was  determined 
that  in  a  worse  case  scenario,  when  for  some  reason  the  dynamic  model  in  EKF  provides 


no  additional  information  to  the  measurements,  then  the  quality  of  attitude  estimation  was 
bound  by  the  errors  of  the  best  sensor.  When  poor  gyros  are  used  this  steady  state  error 

bound  is  the  performance  of  the  interferometer  which  for  a  0.6  m2  array  is  about  0.45 

degrees  1  a  RSS  attitude  error  and  about  0.28  degrees  for  aim2  array  within  the  GPS 
update  range  mentioned  above. 

Varying  gyro  qualities  were  also  analyzed.  For  a  projected  MM  Gyro  with  5°/hr  bias,  8 
hour  time  constant  and  0.25°/hr1/2  ARW,  steady  state  1  a  RSS  attitude  errors  of  about  0.3 

and  0.2  degrees  can  be  expected  for  a  0.6  m2  array  and  1  m2  array,  respectively,  at  a  1/10 
Hz  GPS  update  frequency. 

The  inherent  robustness  to  failure  aspect  of  the  tightly  coupled  integration  scheme  imple¬ 
mented  by  this  system  was  also  analyzed.  Using  both  a  linear  covariance  analysis  and 
actual  space  orbit  simulations,  the  redundancy  of  this  system  was  proven.  Various  sub- 
optimal  cases  were  considered  to  include  loss  of  gyro(s)  and  loss  of  a  GPS  antenna. 

Chapter  6  presented  EKF  results  using  real  gyro  data  obtained  from  a  static  laboratory  test 
with  real  GPS  phase  measurements  collected  on  the  Draper  Laboratory  roof.  Results  con¬ 
firm  that  the  proposed  integration  scheme  is  experimentally  valid  with  current  technology. 

Chapter  7  summarized  the  energy  expenditures  for  the  system.  It  was  concluded  that  an 
ADS  with  1  a  RSS  attitude  errors  of  about  a  quarter  of  a  degree  could  operate  at  about  200 
mWatts.  The  table  below  details  what  a  typical  ADS  of  the  type  needed  to  fill  the  telecom¬ 
munication  niche  costs  in  terms  of  size,  power  and  money.  Note  that  the  following  quotes 
are  rough  estimates  and  were  obtained  from  the  specifications  of  typical  components 
available  in  industry. 
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Sensor 

Mass  (kg) 

Power  (W) 

Cost  (K$) 

INS 

0.75 

11 

60-85 

Sun  Sensor 

0.5 

.5 

40-60 

Horizon  Sensor 

10 

-200 

GPS 

5 

150 

Total 

-6.75 

~27 

-470 

Table  8.1:  Typical  ADS  Components 


The  following  table  summarizes  the  same  type  of  estimate  for  the  proposed  IGPS/MM 
Gyro  ADS. 


Sensor 

Mass  (kg) 

Power  (W) 

Cost  (K$) 

MM  Gyro  (triad) 

0.5 

.15 

50 

IGPS  (1/10  Hz) 

1 

.05 

50 

Total 

~1.5 

-0.2 

-100 

Table  8.2:  Estimates  of  the  IGPS/MM  Gyro  ADS 


Note  that  the  gyro  cost  figure  is  a  rough  estimate  for  a  system  that  is  still  in  its  develop¬ 
mental  states.  From  the  above  tables  it  is  easy  to  see  how  such  a  system  could  find  a  viable 
market  in  the  growing  telecommunications  arena. 

8.2  Suggestions  for  Future  Work 

The  author  can  suggest  several  areas  where  future  work  can  be  concentrated.  Due  to  time 
constraints  for  this  thesis,  a  dynamic  real  data  integration  was  not  possible.  As  stated  in 
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Chapter  6,  MM  Gyro  data  can  be  obtained  from  gyros  mounted  on  turn-tables  pro¬ 
grammed  at  desired  rates.  This  output  can  then  be  integrated  into  the  space  simulation 
used  in  this  analysis.  This  would  serve  as  a  further  step  in  validating  the  system,  and  it 
could  be  used  as  another  way  of  testing  its  sensitivities. 

Another  area  of  study  would  be  research  into  robustness  methods  for  the  model  based  esti¬ 
mation  filter.  Research  in  this  area  has  already  begun  at  Draper,  but  integration  of  robust¬ 
ness  algorithms  with  actual  data  and  a  tightly  coupled  EKF  has  yet  to  be  accomplished. 
Also,  better  dynamic  models  need  to  be  developed.  This  study  only  characterized  gravity 
gradient  effects,  but  other  torques  such  as  internal  spacecraft  torques  and  control  torques 
have  yet  to  be  modeled. 

Finally,  the  GPS  error  models  could  benefit  from  better  modeling  of  phase  center  variation 
and  multipath.  Real  time  identification  of  multipath  has  yet  to  be  successfully  imple¬ 
mented  in  such  a  design.  Whether  the  benefits  obtained  from  such  calibration  justify  the 
cost  (especially  in  a  low  multipath  space  environment)  is  another  question  left  to  be 
answered. 
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Appendix  A 


Derivation  of  Attitude  Rate  Dynamics 

A.l  Introduction 

As  stated  in  Chapter  3,  the  model  based  nature  of  the  Extended  Kalman  filter  used  in  this 
study  is  highly  dependent  on  the  quality  of  modelling  of  the  dynamics.  The  following  will 
give  a  derivation  of  the  dynamics  of  attitude  rate  error  state,  8cq5W  .  Following  the  deriva¬ 
tion  will  be  a  brief  explanation  of  how  this  model  is  incorporated  into  the  framework  of 
the  EKF. 


A.2  Derivation 

The  angular  rate  error,  8cofiA, ,  is  the  difference  between  the  true  rate  and  estimated  rate: 

B  B  B 

^-BN=  -BN  ~  ®BN  (A.l) 

To  derive  the  dynamics  of  the  angular  rate  error  we  first  begin  with  Euler’s  vector  equation 
of  motion: 


r'(T-mBB,x(imBB,))  <a.2) 

Applying  equation  A.  1  we  can  linearize  A.2  to  yield 

^(6  +  5  (Q)h  =  (/  1{r  +  8r-(w  +  6©)g/x(/(w  +  6o))g/)})  (A.3) 

Subtracting  the  estimate  of  A.2  and  ignoring  higher  order  terms  gives 

^(ScoB/)=  I  -  [®s/ x  (A .4) 


Using  equation  3.42  we  can  expand  the  above  to  yield: 


d  /Sj  b 
di^~BN 


)  =  rWf,] 


Bnx 


[a )BBIxI))(b(DB 


BN  +  5®A'/)  +  1 


5  T  (A.5) 


By  remembering  the  relation: 
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(A.6) 


we  know  that: 


d  /  p  B  .  d  .  o  N  v  o  B  B 

^(5ww)  =  ^(5fflw)  -  6©w  x  Sco^  (A.7) 

Furthermore,  we  can  eliminate  the  second  term  above  since  for  a  circular  orbit  co^y  is  con¬ 
stant  leading  A.5  to  be  expressed  as: 

ft(S ISBBN)  =  r1([/mB;]x-[ffl“/Xn)(8ffl®A,  +  5oj“/)-6o!®,x8ffl®„  +  rl5r  (A.8) 

Torque  error  is  comprised  of  contributions  from  the  gravity  gradient  and  other  sources  and 
can  be  expressed  via  the  small  angle  errors  [23]  as 


dTc 

30 


5T  =  ^-gg50  +  — gg50  +  5r 


3TC 

30 


other 


(A.9) 


The  second  term  on  the  right  side  is  zero  since  gravity  gradient  torque  is  a  function  only  of 
angle  displacement  and  not  the  rate.  We  ignore  the  contributions  from  other  torque  sources 
as  well.  In  accordance  with  the  analysis  performed  by  Puri  [20] ,  estimating  torque  errors 
such  as  aerodynamic  or  solar  contributions  yielded  minimal  improvements  in  a  coupled 
interferometric-gyro  attitude  estimator  as  long  as  the  process  noise  on  the  gyro  replace¬ 
ment  state  was  increased  accordingly.  The  remaining  torque  error  comes  entirely  from  the 
gravity  gradient  and  is  given  without  derivation  [23]  as: 


3  Trr 
5  T  =  ^—GG 
30 


50  = 


[6(0 Ni([rx]I-[Ir]x)[rx^BN 


(A.  10) 


where  the  latter  representation  assumes  a  circular  orbit.  We  finally  have  an  expression  for 
the  angular  rate  error  dynamics  from  the  nav  to  body  frame: 
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(A.l  1) 


=  r‘[6ffl2m([fx]/-[/r]x)[^]]84^ 

+ r'ajfliJ,]*  -  taj,  x  /])(SidJw + -  SB®,  x  S(sbbn 

The  linearized  attitude  dynamics  then  fall  out  as  presented  in  Chapter  3.  Following  a  simi¬ 
lar  derivation  the  pure  attitude  rate  dynamics  which  are  used  to  generate  truth  for  the 
Extended  Kalman  filter  follow  the  equation  [20]: 

®BN  =  1  [^-©S/X  (A.12) 

B 

The  term  h(Qm  is  found  by  realizing  that 

—NI  ~  &N1  =  (^BN  '  ®At/)  ~  (CbN  '  =  (C  BN  ~  (A.13) 

and  remembering  that 

Cbn  =  CbnCbn  ~  [1  +  [8©]*]  (A.  14) 

which,  when  applied  to  A.13  yields 

B  B 

S%/  =  [8@]  X'®NI  (A.  15) 
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Appendix  B 


Roof-top  Test  Conditions 

The  conditions  for  the  roof-top  data  collection  of  January  14,  1998,  are  summarized  below.  Time 
is  in  seconds  and  givein  in  GPS  time.  The  differential  phase  measurements  on  this  day  used  the 
GPS  satellite  geometries  indicated  on  the  first  figure.  The  asterisk  shows  when  the  SV  trajectory 
began.  The  second  figure  shows  what  channels  were  used  at  specific  times. 


LOS  Vectors  forOll^gl 


Figure  B.l:  GPS  Satellite  Geometries  During  Data  Collection 
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Satellites  on  Channels  for  011401 
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Figure  B.2:  Channels  for  Data  Take 

The  thicker  line  regions  are  caused  by  circles  being  graphed  instead  of  a  line  and  correspond  to 
times  when  there  was  an  invalid  measurement.  However,  all  data  for  these  regions  is  not  bad  since 
an  invalid  measurement  at  one  second  causes  a  circle  to  be  placed  on  the  graph  which  takes  up 
more  space  than  the  second  for  which  it  corresponds.  The  following  graph  shows  the  total  ADOP 
for  the  geometries  of  this  data  take. 
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deg  /  mm 


SD  KF  ADOP  for  01 1401 


PDOP  for  01 1401 

Figure  B.3:  Total  ADOP  for  Data  Take 
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Appendix  C 


IGPS  Particulars 

C.l  Introduction 

The  following  Appendix  contains  a  brief  derivation  of  Attitude  Dilution  of  Precision 
(ADOP)  and  a  description  of  the  single  difference  (SD)  operator.  More  information  on 
these  two  terms  can  be  found  in  [20,  23], 


C.2  ADOP 

Dilution  of  Precision  (DOP)  is  a  term  that  describes  how  SV  geometry  degrades  accuracy. 
Like  the  more  common  Geometric  Dilution  of  Precision  (GDOP),  ADOP  provides  a  mea¬ 
sure  of  the  theoretical  attitude  solution  accuracy  possible  with  a  given  satellite  geometry. 
Recall  from  equation  3.61  that  the  sensitivity  of  the  differential  phase  to  attitude  error  is: 


3z.  :  B 

-fir  =  2.[p  hest]x.bio 

95q  BN 


DOP  has  the  general  formula  described  by  equation  C.l: 

DOP  =  VtraceK#7#)-1] 

so  for  the  attitude  case, 


(3.62) 


(C.l) 


ADOP  =  ^rnce[(H^69)_1; 


(C.2) 


For  the  ADOP  calculations  in  this  thesis,  however,  ADOP  values  correspond  to  the  follow¬ 
ing  equation: 


ADOP  =  ^  +  <70  + <4 

which  is  effectively  twice  the  value  computed  in  equation  C.2. 


(C.3) 
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C.3  Single  Difference  Operator 

When  a  differential  phase  measurement  is  taken,  a  slave  antenna  phase  is  subtracted  from 
the  master  antenna  phase.  This  procedure  is  represented  by  the  single  difference  operator: 
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GPS  errors,  such  as  multipath,  phase  center  variation,  receiver  noise  and  line  bias,  are 
errors  which  corrupt  the  carrier  phase  before  the  antenna  phase  differencing  takes  place. 
(Baseline  length  error  is  an  exception  since  it  does  not  effect  the  raw  carrier  phase.)  The 
resulting  differential  phase  errors  (for  a  three  baseline  example)  are  represented  by  the  sin¬ 
gle  difference  operator: 
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where  e  is  an  error  vector  on  undifferenced  carrier  phase,  and  Ae  is  the  effect  of  the  error 
on  differential  phase.  The  correlation  between  baselines  is  described  by: 


Aa£  =  E[ Ae  •  A/]  =  E[SDe  ■  eTSDT]  =  SD  •  Ae  •  SDT 


2  1  1 
1  2  1 
1  1  2 


(C.3) 


where  is  the  covariance  of  the  error  on  undifferenced  carrier  phase.  The  differential 
phase  error  covariance,  AAe ,  has  two  times  the  variance  of  the  original  error  and  is  corre- 
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lated  across  baselines  [20]. 
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